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1.  Introduction 


This  work  was  a  collaborative  effort  between  Duke  University  and  SAIC  (formerly  AETC 
Incorporated).  The  final  report  consists  of  two  volumes.  This  first  volume  gives  a  complete 
description  of  the  overall  effort  as  it  evolved  and  culminated  in  a  successful  demonstration  of 
IMU-based  positioning  of  a  handheld  EM61  sensor  for  data  collection  to  support  target 
characterization  of  buried  UXO.  Volume  2  is  a  complete  report  of  the  supporting  activities  at 
Duke  University. 

1.1.  Background 

Past  studies  have  clearly  demonstrated  that  commercial  UXO  sensors  (specifically,  the  Geonics 
EM61  and  EM61-HH,  and  the  Geophex  GEM-3)  can  be  used  effectively  to  discriminate  between 
buried  UXO  and  clutter  if,  and  only  if,  the  sensor  position  and  orientation  are  precisely  known 
while  the  data  are  collected  ([  1  ]— [4]).  The  problem  is  that  it  is  not  always  easy  to  get  sensor 
location  information  at  the  levels  of  accuracy  needed  for  reliable  discrimination  and 
classification,  i.e.  within  about  a  cm.  A  simple  solution  to  the  sensor  positioning  problem  is  to 
collect  data  on  a  fixed  grid  over  the  target  using  a  template  laid  on  the  ground  (see  Figure  1). 
Although  inexpensive  and  robust,  this  approach  leaves  much  to  be  desired.  A  positioning  system 
that  is  integral  to  the  sensor  and  allows  the  operator  to  sweep  the  sensor  around  above  the  target 
is  much  preferred.  This  would  permit  use  in  rougher  terrain  and  confined  locations.  It  would  also 
facilitate  collection  of  a  much  higher  density  of  data  over  the  target.  The  grid  approach  provides 
a  limited  number  of  samples  (e.g.  36  for  a  6x6  grid)  and  requires  signal  to  noise  ratios  greater 
than  25  or  30dB  in  order  to  get  good  target  classification.  Higher  data  density  translates  directly 
to  improved  signal  to  noise  ratio  at  the  processor  output,  and  hence  improved  performance  for 


Figure  1:  Collecting  EM61-HH  data  on  a  grid  over  a  suspected  buried  UXO  item. 
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deeper  targets.  It  is  impractical  to  try  to  improve  the  data  density  with  a  grid  template  approach. 
However,  with  a  freely  swept  system  operating  at  10Hz,  the  data  density  can  be  increased  by  a 
factor  of  50-80  during  the  time  required  to  carefully  sample  the  grid  using  the  current  approach 
(3-5  minutes). 

1.2.  Objective 

The  goal  of  this  project  was  to  develop  an  inexpensive,  robust  way  to  accurately  determine  the 
trajectory  of  a  handheld  UXO  sensor  as  it  is  swept  about  above  a  suspected  buried  UXO  item.  In 
principle,  the  sensor  position  and  orientation  can  be  tracked  using  an  inertial  measurement  unit 
(IMU)  that  measures  accelerations  and  angular  rates  along  three  orthogonal  axes.  Whether  a 
compact,  inexpensive,  rugged  IMU  based  on  micro-electromechanical  systems  (MEMS) 
technology  is  sufficient  to  do  the  job  was  a  main  focus  of  this  project.  Another  point  of  focus 
was  the  development  of  processing  procedures  that  combine  the  trajectory  with  the  sensor  output 
stream  to  characterize  the  target  for  discrimination  between  UXO  and  clutter 


1.3.  Technical  Approach 

The  approach  taken  early  stages  of  this  project  is  summarized  in  the  block  diagram  below. 


The  first  step  was  to  characterize  and  develop  parametric  descriptions  for  the  sensor  motion  and 
IMU  errors,  followed  by  the  development  of  adaptive  processing  and  filtering  schemes  to  reduce 
the  IMU  errors.  The  idea  here  being  to  use  the  IMU  errors  to  drive  the  estimation  of  the  state  of 
the  sensor  motion.  Within  this  scheme,  external  aids  to  ensure  that  the  IMU  and  position  errors 
remain  strongly  correlated  and  converge  to  zero  -  such  as  the  UXO  sensor  output  or  other 
independent  references  -  were  to  be  explored. 
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Initially,  the  project  proceeded  according  to  plan.  Tests  were  conducted  with  the  IMU  and  the 
EM61-HH  in  controlled  laboratory-type  environments.  Significantly,  this  early  work  focused  on 
constrained  two-dimensional  motion,  and  it  was  found  that  with  rather  straightforward 
processing,  IMU-based  positioning  could  be  obtained  with  the  precision  required  to  support 
accurate  dipole  inversion  of  EM61-HH  data.  The  IMU  unit  used  in  this  early  work  was  an 
inexpensive  Crossbow  IMU400CC-100.  The  data  collections  are  described  in  §2,  and  the  basic 
processing  approach  for  integrating  the  IMU  data  is  described  in  §3.  In  order  to  get  accurate 
dipole  inversion  of  the  EM61-HH  data,  we  had  to  incorporate  a  forward  model  that  included  the 
dynamic  response  characteristics  of  the  EM61-HH,  as  discussed  in  §4. 

During  this  period,  the  Duke  University  research  effort  focused  on  several  signal  processing 
approaches  for  handheld  UXO  sensor  improvements  that  could  be  incorporated  in  the  overall 
processing  flow,  to  begin  they  explored  system  identification  as  a  means  to  directly  model  the 
errors  in  the  measured  quantities  (acceleration  and  angular  rate).  In  this  application,  system 
identification  determines  the  transfer  function  of  the  unknown  filter  from  known  input  and 
output  data,  thus  modeling  the  relationship  between  the  observed  IMU-derived  data  (output)  and 
the  ground  truth  (input).  Because  the  IMU  accelerometers  and  angular  rate  sensors  are  most 
accurate  during  different  types  of  motion  and  have  different  limitations,  adaptive  error  mitigation 
algorithms  employing  feedback  loops  to  use  the  more  accurate  sensor  (accelerometer  or  angular 
rate)  to  mitigate  errors  in  the  less  accurate  sensor  (accelerometer  or  angular  rate)  were  then 
developed  and  tested.  This  effort  resulted  in  iterative  integration  algorithms  for  velocity 
stabilization  and  position  calculation  which  were  found  to  perform  well  for  two-dimensionally 
constrained  sweeping  motions.  A  zero  velocity  update  (ZUPT)  algorithm  was  introduced  for 
improved  calculations  of  the  velocities.  The  yaw  angular  rate  was  used  to  determine  the  zero- 
velocity  points  for  nonlinear  motion,  and  the  acceleration  and  its  integration  quantity  to  find  the 
zero-velocity  points  for  linear  motions.  Results  obtained  for  the  two-dimensionally  constrained 
motion  tests  indicated  that  this  approach  is  capable  of  accurately  positioning  the  handheld  UXO 
sensor. 

We  used  video  capture  for  ground  truth  with  the  two-dimensionally  constrained  tests.  That 
approach  was  not  viable  for  unconstrained  three  dimensional  motions.  The  Arcsecond  laser 
positioning  system  was  evaluated  and  found  to  be  acceptable  for  ground  truth  for  general  three- 
dimensional  motions,  as  discussed  in  §5.  The  tests  with  three  dimensional  motion  are  described 
and  discussed  in  §6.  The  algorithms  that  worked  well  with  the  two-dimensionally  constrained 
motion  failed  to  produce  accurate  position  information  with  the  Crossbow  IMU400CC-100  for 
three-dimensional  motion.  This  was  borne  out  by  the  Duke  results.  When  applied  to 
unconstrained  three  dimensional  motions,  estimates  of  the  system  parameters  failed  to  converge, 
and  plots  of  the  measured  accelerations  versus  the  true  accelerations  showed  no  clear 
relationship  between  the  measured  and  true  accelerations.  Overall,  processing  of  unconstrained 
three-dimensional  data  produced  results  which  do  not  meet  the  desired  sensor  positioning 
performance. 

We  concluded  that  the  reason  why  we  had  good  results  in  detennining  sensor  location  from  the 
Crossbow  IMU  data  for  2D  motion  but  discouraging  results  for  3D  motion  was  because  of 
unacceptable  levels  of  noise  and  drift  in  the  Crossbow  angle  rate  measurements.  Accurate 
determination  of  the  IMU's  orientation  is  vital  to  getting  accurate  position  information  because 
even  a  very  small  error  in  angle  can  introduce  erroneous  gravity  contributions  to  the  X  and  Y 
accelerations,  which  in  turn  integrate  up  to  large  position  errors.  A  higher  quality  Honeywell 
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HG1900  IMU  was  integrated  with  the  EM61-HH  by  ENSCO  under  ESTCP  project  MM-0604.  It 
is  referred  to  as  the  Small-Area  Inertial  Navigation  Tracking  or  SAINT  system.  We  collaborated 
with  MM-0604  to  evaluate  whether  or  not  the  SAINT  could  provide  sufficiently  accurate  sensor 
position  infonnation  to  support  reliable  UXO/clutter  discrimination.  As  discussed  in  §7,  we  were 
able  to  get  discrimination  quality  geolocated  data  using  the  SAINT  system  at  the  Aberdeen 
Standardized  UXO  Test  Site. 

Finally,  in  §8  we  summarize  work  demonstrating  the  feasibility  of  including  sensor  data  in  the 
track  estimation  process.  The  approach  is  based  on  the  principle  that  the  best  model/data  fit 
occurs  when  the  sensor  trajectory  is  correct.  Slight  errors  in  the  trajectory  degrade  the  fit  quality. 
The  idea  is  then  to  assume  that  the  trajectory  from  the  IMU  is  in  error,  and  allow  slight 
modifications  of  the  trajectory  in  the  EMI  inversion  process,  iterating  until  the  best  fit  to  the  data 
is  achieved.  Trajectory  errors  typically  arise  because  of  errors  in  estimating  the  IMU  biases.  By 
including  the  bias  parameters  in  our  EM  inversion  procedure,  we  were  able  to  simultaneously 
determine  the  best  track  estimate  along  with  the  target  parameters  (location  and  principal  axis 
polarizabilities).  Simulations  using  this  procedure  showed  that  we  can  reduce  5  cm  rms  track 
errors  to  <5  mm  nns. 
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2.  Controlled  Two-Dimensional  Motions 


The  program  plan  called  for  a  sequential  approach  using  increasingly  complex  sensor  motions: 

a)  Constrained  Path  -  Tightly  constrained  motion,  with  the  sensor  level  and  following  a 
simple  zigzag  path  above  the  target; 

b)  Controlled  Sweep  -  A  natural  sweeping  motion  back  and  forth  over  the  target  where  the 
orientation  of  the  sensor  head  varies  in  a  natural  fashion;  and 

c)  Unconstrained  Motion  -  Arbitrary  sensor  path  above  the  target. 

This  section  describes  the  data  collection  process  for  controlled  two-dimensional  motions, 
comprising  a)  and  b)  above. 

2.1.  Measurement  Setup 

2.1.1.  Apparatus  Description 

The  apparatus  used  in  the  data  collection  process  is  shown  in  Figure  2  below.  A  Crossbow 
IMU400CC-100  and  a  Geonics  EM61-HH  were  secured  at  known  fixed  positions  on  a  rigid 
wooden  board.  Collectively,  the  board  and  sensors  are  referred  to  throughout  this  report  as  the 
body.  It  is  this  body  that  we  are  interested  in  carefully  tracking  -  initially  through  various  tightly 
constrained  maneuvers,  and  later  by  gradual  easing  of  the  constraints,  as  the  above  list  details. 


Figure  2:  Two  views  of  the  body  consisting  of  a  Crossbow  IMU400  (with  its  sensitive  axes  expanded  on  the 
left)  and  a  Geonics  EM61-HH. 


The  Crossbow  IMU400CC-100  is  a  solid-state  six-degree-of-freedom  (6DOF)  Inertial  System 
designed  for  general  measurement  of  the  acceleration  and  angular  rate  vectors  in  dynamic 
environments.  Its  physical  dimensions  are  7.62  x  9.53  x  8.13  cm  and  it  weighs  less  than  640  g. 
The  three  angular  rate  sensors  are  bulk  micro-machined  vibratory  MEMS  sensors  that  utilize  the 
Coriolis  force  to  measure  the  angular  rates  along  each  of  the  sensitive  axes  labeled  in  Figure  2. 
The  three  MEMS  accelerometers  are  surface  micro-machined  silicon  devices  that  use  differential 
capacitance  to  independently  sense  the  acceleration  components  along  these  same  axes.  When 
operating  in  the  continuous  update  mode,  the  sampling  rate  of  the  IMU  is  just  below  135  Hz. 
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The  Geonics  EM61-HH  is  a  handheld  time  domain  or  pulsed  electromagnetic  induction  (EMI) 
sensor  with  the  transmit  and  receive  coils  offset  from  each  other.  The  transmitted  primary  field 
induces  eddy  currents  in  the  target,  which  in  turn  radiate  a  secondary  field  measured  by  the 
receive  coil.  The  EM61-HH  measures  the  secondary  field  at  time  delays  of  147,  263,  414  and 
613  ps  after  the  primary  field  shutoff  at  a  sampling  rate  of  15  Hz. 


Figure  3:  The  apparatus  as  viewed  from  an  overhead  video  camera.  The  point  of  attachment  also  served  as 
a  pivot  point. 


As  is  apparent  from  Figure  3  above,  the  body  was  laid  on  another  rigid  wooden  board  and 
attached  to  this  at  a  point.  The  point  served  as  a  pivot  point,  about  which  the  body  was 
constrained  to  rotate  in  a  purely  planar  trajectory.  In  addition,  the  underlying  board  was  set  on  a 
linear  rail  system  offering  an  additional  translational  component  to  the  motion  of  the  body, 
thereby  allowing  the  body  to  move  in  a  variety  of  planar  patterns,  including  the  simple  zigzag 
paths  alluded  to  in  the  list  above. 

A  video  camera  was  installed  overhead  the  apparatus  to  provide  ground  truth.  In  fact,  Figure  3  is 
a  single  frame  extracted  from  the  video  data  for  one  of  several  measurement  scenarios 
undertaken.  As  a  calibration  aid  for  the  video,  a  peg  board  with  holes  spaced  1”  apart  in  a  regular 
2D  grid  pattern  was  used  as  the  background.  A  smaller  piece  of  the  same  board  with  4  clearly 
labeled  stickers  was  also  attached  to  the  EM61-HH  sensor  head  (see  Figure  4).  These  stickers 
were  essential  in  tracking  the  body  movement  frame-by-frame,  thereby  allowing  the  extraction 
of  ground  truth  positioning  information  from  the  video  data.  The  extraction  process  was  done 
manually,  frame  by  frame,  using  a  set  of  IDL  routines  to  display  the  images  and  record  locations 
of  the  cursor  when  centered  over  the  dots. 
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Figure  4:  Stickers  placed  on  the  EM61-HH  sensor  head  for  ground  truth  positioning  purposes. 
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2.1.2.  Performed  Scenarios 


Several  measurement  scenarios  were  conducted.  Figure  5  -  Figure  9  below  show  the  variety  in 
these  scenarios,  spanning  from  a  simple  linear  trajectory  (Figure  5)  to  the  more  complicated 
zigzag  sweeping  trajectory  (Figures  8  or  9).  Since  the  mechanism  of  extracting  ground  truth 
positioning  information  from  the  video  data  via  the  stickers  proved  to  be  labor  intensive,  only  a 
limited  portion  of  the  video  data  for  each  measurement  scenario  undertaken  was  subjected  to  this 
process.  In  fact,  Figure  5  -  Figure  9  represent  all  the  ground  truth  positioning  information  that 
was  extracted.  Fortunately,  this  was  more  than  enough  data  for  the  sensor  motion 
characterization  process. 


Figure  5:  Data  Set  A  -  Back  &  Forth  Linear  Trajectory  over  a  2”  A1  sphere. 


Figure  6:  Data  Set  B  -  Back  &  Forth  Arced  Trajectory  over  a  2”  A1  sphere. 
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Figure  7:  Data  Sets  Cl  &  C2  -  Two  Back  &  Forth  Sweep  and  Push  Trajectories  over  a  2”  A1  sphere. 


Figure  8:  Data  Set  D  -  Back  &  Forth  Swing  Trajectory  over  a  2”  A1  sphere. 
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Figure  9:  Data  Set  E  -  Back  &  Forth  Swing  Trajectory  over  a  40mm  Projectile. 


2.2.  Definitions  and  Notation 

This  section  introduces  the  notation  used  throughout  this  report  and  defines  some  key  elements 
needed  in  comparing  the  IMU-derived  results  to  the  video  ground  truth. 

Following  standard  practice,  all  vectors  and  matrices  are  represented  as  bold  lowercase  and 
uppercase  letters,  respectively.  Each  vector  contains  a  superscript  letter  specifying  the  coordinate 
frame  in  which  the  vector  is  defined,  as  well  as  an  optional  subscript  identifying  an  attribute  of 

the  vector  (e.g.  rsv  represents  the  position  vector  of  sticker  4  as  defined  by  the  video  coordinate 

frame,  while  a*M£/  represents  the  IMU  acceleration  vector  as  defined  by  the  body  coordinate 
frame). 

In  general,  to  transfonn  a  vector  from  one  coordinate  frame  to  another  requires  a  combination  of 
a  translation  and  a  rotation.  The  rotation  can  be  completely  defined  by  an  orthonormal  matrix 
representing  an  ordered  sequence  of  three  plane  rotations  involving  the  Euler  angles  -  yaw,  pitch 
and  roll.  An  explicit  representation  of  the  rotation  matrix  R6o_>6(0  needed  to  transform  a  vector 

defined  by  the  initial  body  frame  at  time  t  =  to  to  the  body-frame  at  time  t  is  developed  in 
Appendix  A.  Similar  representations  hold  for  the  rotation  matrix  Rv^i(i)  needed  to  transform  a 
vector  defined  by  the  video  frame  to  the  body-frame  at  time  t  and  the  rotation  matrix  Re^.fe(t) 
needed  to  transform  a  vector  defined  by  the  earth  frame  to  the  body-frame  at  time  t,  the 
difference  being  in  the  actual  angle  values. 

Referring  to  Figure  10  below,  it  is  clear  that  if  the  IMU  represents  the  origin  of  the  body 
coordinate  frame  Ob  with  the  IMU  sensitive  axes  coinciding  with  the  body  axes,  then  in  the 
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video  coordinate  frame  the  position  vector  of  the  z'-th  sticker  at  time  t  -  i.e.  r‘(7)  -  can  be 
expressed  in  terms  of  the  IMU  position  vector  at  time  t  as 

(1)  r’(0  =  r^(0  +  R^(0<  ^(0  +  R^R^(f)< 

Here,  d  *  represents  the  position  vector  of  the  z'-th  sticker  as  defined  by  the  body  coordinate 
frame,  and  RA_>V(0  =  R.  >VR6.  >hJt)  represents  the  rotation  matrix  that  transforms  d*  from  the 
body  frame  to  the  video  frame  in  terms  of  a  rotation  matrix  that  first  transforms  d/'  from  the 
body  frame  at  time  t  to  the  initial  body  frame  at  time  t=to,  followed  by  a  time  invariant  rotation 
matrix  that  transforms  d!r  from  the  initial  body  frame  to  the  video  frame. 


Figure  10:  Schematic  diagram  showing  the  relevant  vectors  and  coordinate  frames. 


It  should  also  be  clear  from  Figure  10  that 

rjMU  (0  =  Cu  (t  =  K )  +  R,^  (0 
so  that  (1)  can  be  rewritten  as 
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(2)  !•;;  (0  =  x]MU  (t  =  t0)  +  R^v  ( r^j  (0  +  (0< ) 

Finally,  noting  that 


c,(^=o=r;l(^=o-R^v<: 

we  have,  upon  substituting  into  (2), 

(3)  r'(0  =  r’(<  =  O  +  R^bH„W-dS  +R,^(0d‘ ) 

Equation  (3)  is  the  necessary  expression  whereby  the  IMU-derived  results  can  be  compared  to 
the  video  ground  truth. 
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3.  Data  Analysis  and  Signal  Processing 


3.1.  Inertial  Measurement  Unit 
3.1.1.  Attitude  Determination 


The  attitude  of  the  IMU  (and  thus  the  body)  is  detennined  from  the  measured  angular  rates  cox , 
coy  and  coz .  Based  on  the  Euler  angle  rotations,  t//,  #  and  (f)  (yaw,  pitch  and  roll,  respectively) 
described  in  Appendix  A,  we  have  upon  referring  to  Figure  Al, 


d6 

y  ■  and 


where  x",  y"  and  z'  are  the  unit  vectors  about  which  the  respective  rotations  take  place. 
Transforming  these  vectors  to  the  body-frame  at  time  t,  we  have  from  (A3) 
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Thus,  we  can  write 


“1 

0 

-  sin# 

= 

0 

COS  (j) 

cos#  sin  (j) 

0 

-  sin^ 

COS#  COS  (j) 

from  which  we  get,  by  inversion, 
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For  a  level  body  constrained  to  move  in  a  plane,  9  =  </> 


0  and  so  (4)  results  in  cox 


d(j) 

dt 
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dt 
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This  means  that  the  only  non-zero  angular  rate  component  is  along 


the  body  z-axis,  as  expected,  and  the  yaw  angle  can  be  determined  by  a  straightforward 
integration  of  coz .  When  d±  0  and/or  $  *  0  ,  however,  the  angular  rates  no  longer  decouple  and 
the  angles  now  need  to  be  determined  via  (4). 


Although  every  effort  in  the  measurement  process  was  made  to  constrain  the  body  to  move  in  a 
plane,  deviations  from  this  ideal  inevitably  creeps  in.  As  an  example,  the  IMU  angular  rates  for 
Data  Set  C2  are  shown  below  in  Figure  1 1.  It  clear  from  the  figure  that  both  the  x  and  y  angular 
rate  components  are  non-zero.  This  suggests  that  (4)  must  be  used  if  the  goal  is  to  accurately 
characterize  the  motion  of  the  body. 

All  three  angular  rates  contain  biases  and  trends.  The  red  lines  in  Figure  1 1  represent  linear  fits 
to  the  t  <  -7s  and  t  >  50s  data,  when  the  body  was  known  to  be  stationary  and  in  the  same  general 
position.  The  respective  bias  and  drift  values  of  the  three  angular  rate  components  are  listed  in 
Table  1  for  this  data  set,  as  well  as  several  other  data  sets. 


Data 

Set 

(Ox  bias 
(%) 

(Ox  drift 

(xl0'57s2) 

(ov  bias 
(7s) 

(ov  drift 

(xl0‘57s2) 

(oz  bias 
(7s) 

(Oz  drift 

(xl0'57s2) 

Cl 

0.1538655 

13.43 

0.0709080 

-5.94 

0.0900197 

12.51 

C2 

0.1605218 

3.25 

0.0690165 

15.67 

0.0952184 

-6.87 

D 

0.1657144 

-19.74 

0.0849374 

14.80 

0.0911045 

20.19 

E 

0.1769275 

18.67 

0.0674774 

-3.56 

0.0855456 

23.50 

Tablel:  Biases  and  linear  drift  rates  for  the  angular  rates  of  four  separate  data  sets. 


The  angles  calculated  using  (4),  where  the  angular  rates  have  had  their  linear  trends  removed,  are 
shown  in  Figure  12  for  Data  Set  C2.  The  angles  in  Figure  12  have  also  had  a  linear  trend  (again 
obtained  from  the  t  <  -7s  and  t  >  50s  data)  removed  from  them.  The  respective  bias  and  drift 
values  of  the  three  angles  are  listed  in  Table  2  below  for  the  same  data  sets  as  in  Table  1.  The 
blank  cells  in  the  table  represent  instances  where  a  linear  trend  based  on  before  and  after 
stationary  data  is  either  inapplicable  (i.e.  the  before  and  after  data  may  represent  clearly  different 
positions),  or  is  inadequate  to  describe  the  drift. 


Data 

Set 

<f>  bias 
(xlO'3  °) 

(f)  drift 

(xl0'37s) 

9  bias 
(xlO'3  °) 

9  drift 

(x10'37s) 

y/  bias 
(xlO'3  °) 

y/  drift 

(xl0'37s) 

Cl 

13.4 

1.0 

17.0 

3.2 

80.8 

7.3 

C2 

-18.6 

-2.1 

-12.8 

-2.2 

-10.5 

-11.6 

D 

E 

5.8 

-5.3 

Table  2:  Biases  and  linear  drift  rates  for  the  determined  angles  for  the  same  data  sets  as  in  Table  1. 
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Figure  12:  Angles  determined  via  the  IMU  angular  rates  for  Data  Set  C2.  The  over-plotted  red  curve 
represents  the  yaw  angle  determined  from  the  video  data. 
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The  red  yaw  curve  over-plotted  in  Figure  12  represents  the  yaw  angle  detennined  by  using  the 
video  ground  truth  via 


where  |LSJ  is  the  magnitude  of  the  lever  arm  vector  that  points  from  the  pivot  point  to  the  7-th 

sticker  as  shown  in  Figure  10;  xoffset  is  the  x-distance  from  the  initial  pivot  point  to  the  origin  of 
the  video  coordinate  frame;  and  i//offset  is  the  angle  between  the  axis  of  linear  motion  (i.e.  the  path 
the  pivot  point  takes  -  see  Figure  10)  and  the  video  y-axis. 

In  the  determination  of  i//  (?) ,  the  following  assumptions  were  made: 


•  L_s.  |  =  0.907  m; 

•  =  0-032  m; 

•  V ob«  =  1.11°;  and 

•  the  video  and  the  IMU  are  perfectly  synchronized 

If  the  same  above  assumptions  are  made  for  the  other  data  sets,  good  general  agreement  between 
the  IMU  determined  yaw  and  the  video  determined  yaw  seems  to  hold  (see  Figure  13).  The  most 
obvious  mismatch  appears  to  be  in  the  trends  of  the  two  curves,  as  well  as  the  relative  timing 
between  the  video  and  the  IMU.  The  trend  mismatches,  also  apparent  in  the  yaw  of  Figure  12, 
are  more  than  likely  due  to  the  imprecise  representations  of  the  IMU  angular  rate  drifts  as  linear 
drifts.  Fortunately,  as  we  shall  see  in  the  next  section,  a  possibility  of  refining  this  representation 
suggests  itself  based  on  enforcing  consistency  between  what  the  angular  rates  tell  us,  and  what 
the  accelerations  tell  us,  about  the  physics  of  the  particular  situation.  The  relative  timing 
mismatches,  on  the  other  hand,  are  greatly  minimized  by  simply  remapping  the  IMU  time  stamps 
as 

(6)  t IMU  — ^  (l  ^t  ^  IMU  t offset  ) 

where  8,  is  a  scale  factor  and  is  fixed  at  a  value  of  0.005  and  t  „  ,  is  the  relative  time  offset 

t  ojjset 

between  the  IMU  and  the  video,  and  is  given  in  multiples  of  sample  time  intervals. 

It  should  be  noted  that  a  linear  trend  was  removed  from  the  yaw  of  Data  Set  Cl  (refer  to  Table  2 
for  the  coefficient  values),  whereas  no  trends  were  removed  from  the  yaws  of  Data  Sets  D  and  E. 
The  latter  was  simply  due  to  the  fact  that  the  choice  of  initial  and  final  yaw  values  were 
different,  making  the  removal  of  a  linear  trend  based  on  the  before  and  after  stationary  data 
inapplicable.  An  effort  was  consciously  made,  wherever  possible,  to  choose  the  IMU  data 
segments  so  that  trend  removal  was  viable.  However,  this  criterion  was  usually  overridden  when 
given  the  option  of  choosing  longer  and  quieter  stationary  segments  of  IMU  data  with  different 
initial  and  final  values. 
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Figure  13:  The  computed  yaw  angles  for  other  IMU  data  sets.  The  over-plotted  red  curves  represent  the 
yaw  angle  determined  from  the  video  data. 
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3.1.2.  Velocity  and  Position  Determination 


The  velocity  and  position  vectors  of  the  IMU  (and  thus  the  body)  are  determined  from  the 
measured  accelerations  ax ,  a  y  and  a: .  Since  the  pivot  point  is  constrained  to  move  along  the 

axis  of  linear  motion,  as  shown  back  in  Figure  10,  it  is  useful  to  explicitly  consider  this  point  in 
the  analysis.  In  the  initial  body  frame,  we  can  write 

(7)  r,i„(()  =  r>(«)-R,„s.(<)(tls.  -lJ‘  =r;-(«)-R, „,,(<)< 


where  r  *°(t)  represents  the  position  vector  of  the  pivot  point  at  time  t  expressed  in  the  initial 
body  frame;  and  (d  s,  -  L  s.  ] b  =  d is  the  vector  that  points  from  the  IMU  to  the  pivot  point  and  is 
independent  of  i.  By  definition,  r^,(.  (t  =  t0)  =  0  . 


Differentiating  (7)  gives  the  following  result  for  the  velocity  of  the  IMU  at  time  t  expressed  in 
the  initial  body  frame  (notation  adopted  from  [5]): 


(8) 


Vp(t) 
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b - ^br. 


(0«v(0d 
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where 

(9)  -  (rov(0x) 


0  -6) At)  CO y  (t) 

CO  At)  0  -00  At) 

-CO  At)  CO  At)  o 


is  a  skew-symmetric  matrix  that  defines  the  cross  product  vector  operation  and  o >^(0  reads  as 

the  rate  of  angular  rotation  vector  of  the  body  frame  at  time  t  relative  to  the  initial  body  frame  at 
time  t  =  t0,  expressed  in  the  body-frame  at  time  t.  Note  that  coit) ,  co  jt)  and  oo2(t)  are  simply  the 

angular  rate  components  as  measured  via  the  IMU  at  time  t,  relative  to  the  respective 
components  measured  at  time  t  =  to, 

Differentiating  (8)  gives  the  following  result  for  the  acceleration  of  the  IMU  at  time  t  expressed 
in  the  initial  body  frame: 


(10) 


4"(0-iW0 


Kb(0  Kb(t) 


dKbit) 

dt 
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where  the  first  tenn  on  the  right  hand  side  represents  the  linear  acceleration  of  the  pivot  point 
and  is  directed  along  the  axis  of  linear  motion  of  Figure  10;  the  second  tenn  represents  the 
centripetal  acceleration  and  is  directed  towards  the  pivot  point;  and  the  third  tenn  is  due  to  the 
angular  acceleration  and  points  in  the  direction  of  rotation,  orthogonal  to  the  centripetal  vector. 
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Since  &IMU(t)  in  (10)  does  not  include  the  influence  of  the  earth  frame  gravity  vector  g  ,  but  the 
IMU  accelerometer  components  do,  we  write 


(11)  a/W  (0  =  (0  [ a  mu  (0  -  Re^,  (0  g] 


where  &bIMu  ( t )  denotes  the  acceleration  vector  at  time  t  as  measured  by  the  IMU. 
Finally,  equating  (10)  and  (1 1)  gives 


(12) 
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IMU  ' 


(t)  =  Rb.b(t)ab;(t)- 


«L(0nL(0 


dtoyit) 

dt 


di+R^(0g 


Equation  (12)  explicitly  represents  the  IMU  accelerometer  values  in  terms  of  the  four  possible 
contributions  to  the  acceleration:  linear,  centripetal,  angular  and  gravitational. 

At  this  point,  it  is  important  to  reflect  on  the  significance  of  equation  (12).  Ignoring  the  linear 
contribution  to  the  acceleration  for  the  time  being,  it  should  first  be  noted  that  the  centripetal  and 
angular  contributions  are  completely  detennined  via  the  IMU  angular  rates  and  a  known  time- 
invariant  vector  (i.e.  d*  ).  It  should  further  be  noted  that  the  gravitational  contribution  is  also 
determined  via  the  IMU  angular  rates,  albeit  indirectly  through  the  body  attitude  (i.e.  angular) 
computations.  However,  here  the  use  of  some  initial  stationary  accelerometer  data  is  also 
necessary  in  order  to  initialize  the  relative  alignment  of  the  earth  and  body  frames.  All  in  all,  the 
essence  of  what  emerges  from  examining  (12)  is  that  independently  measured  quantities  can  be 
related  through  the  physics  of  the  particular  situation,  leading  to  opportunities  for  comparison 
and  correction  schemes. 

In  order  to  elaborate,  we  focus  on  three  specific  well-defined  measurement  situations:  the  pure 
rotational  motion  of  Data  Set  B  (Figure  6);  the  rotational  and  linear  motion  of  Data  Set  C2, 
where  the  motions  are  uncoupled  (Figure  7);  and  the  rotational  and  linear  motion  of  Data  Set  E, 
where  the  motions  are  coupled  (Figure  9). 

3. 1.2.1.  Pure  Rotational  Motion 

Data  Set  B  represents  a  special  case  where  the  linear  contribution  in  ( 12)  is  non-existent.  The  x,  y 
and  z  IMU  accelerometer  values  for  this  data  set  are  shown  in  Figure  14  -  Figure  16.  Within 
each  of  the  figures,  the  right  hand  side  of  (12)  -  as  computed  using  the  IMU  angular  rates  -  is 
shown  in  the  top  panel  in  red,  while  the  separate  centripetal,  angular  and  gravitational 
contributions  are  displayed  in  the  lower  panel.  It  is  clear  from  the  figures  that  there  is  reasonable 
agreement  between  what  the  accelerometers  tell  us  and  what  the  angular  rates  tell  us,  but  that 
there  is  room  for  improvement.  This  is  especially  evident  when  looking  at  the  gravitational 
contribution  in  Figure  14  and  noting  the  bias  differences  relative  to  the  accelerometer  data  during 
stationary  periods  -  i.e.  when  all  the  other  contributions  are  zero.  Also  evident  when  looking  at 
the  angular  acceleration  contributions  in  Figure  15  and  Figure  16,  is  the  fact  that  the 
differentiation  operation  injects  unwanted  noise.  This  problem  can  largely  be  overcome  by  low 
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Figure  14:  IMU  x-component  accelerometer  values  for  Data  Set  B.  The  red  curve  in  the  top  panel 
represents  the  x-component  of  the  right  hand  side  of  equation  (12),  as  computed  via  the  IMU  angular  rates. 
The  colored  curves  in  the  bottom  panel  show  the  respective  contributions  to  the  x-component  of  the 
acceleration. 
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Figure  15:  IMU  y-component  accelerometer  values  for  Data  Set  B.  The  red  curve  in  the  top  panel 
represents  the  y-component  of  the  right  hand  side  of  equation  (12),  as  computed  via  the  IMU  angular  rates. 
The  colored  curves  in  the  bottom  panel  show  the  respective  contributions  to  the  y-component  of  the 
acceleration. 
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Figure  16:  IMU  z-component  accelerometer  values  for  Data  Set  B.  The  red  curve  in  the  top  panel 
represents  the  z-component  of  the  right  hand  side  of  equation  (12),  as  computed  via  the  IMU  angular  rates. 
The  colored  curves  in  the  bottom  panel  show  the  respective  contributions  to  the  z-component  of  the 
acceleration. 
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more  natural  solution  would  be  to  ask  whether  it  is  possible  to  avoid  using  the  differentiation 
operation  altogether  and  still  have  a  viable  error  correcting  scheme. 

To  determine  the  velocity  and  position  of  the  body  for  Data  Set  B,  we  begin  by  subtracting  out 
the  gravitational  contribution  from  the  three  IMU  accelerometer  components.  Figure  17  shows 
the  result  of  integrating  the  gravitationally-compensated  x  and  y  accelerometer  values  to  get  the 
respective  velocity  components  in  the  initial  body  frame  (black  curves).  Also  shown  on  the  same 
plots  (green  curves)  are  the  velocity  components  computed  independently  using  the  IMU  angular 
rates  via  (8).  It  is  obvious,  upon  comparison,  that  the  use  of  the  IMU  angular  rates  produce  much 
more  stable  results.  This  observation  is  compounded  when  looking  at  the  detennined  IMU 
positions  as  shown  in  Figure  18.  In  this  case,  the  black  curves  represent  the  result  of  integrating 
the  black  curves  of  Figure  17,  whereas  the  green  curves  are  computed  directly  from  (7)  using  the 
IMU  angular  rates. 

To  demonstrate  just  how  stable  the  IMU  angular  rate  results  are,  we  compare  the  position 
determination  using  the  IMU  angular  rates  via  (7)  with  the  ground  truth  positions  extracted  from 
the  video  data.  The  comparison  is  given  in  Figure  19  in  the  video  frame  (i.e.  the  coordinate 
system  of  Figure  6)  for  sticker  4,  with  the  ground  truth  data  over-plotted  in  red.  Figure  20  shows 
details  of  the  mismatch  for  both  the  x  and  y  positions.  Although  biases  remain  that  are  likely 
attributable  to  video  distortion  errors,  the  good  news  seems  to  be  that  over  the  length  of  30 
seconds  the  linearly-compensated  measured  angular  rates  appear  to  be  very  stable. 

3. 1.2.2.  Uncoupled  Rotational  and  Linear  Motion 

Data  Set  C2  represents  the  case  where  the  motion  of  the  body  is  decoupled  into  a  rotational 
segment  followed  by  a  linear  segment,  in  alternating  fashion.  The  x,  y  and  2  IMU  accelerometer 
values  for  this  data  set  are  shown  in  Figure  21  and  Figure  22.  In  Figure  21,  the  linear 
contributions  -  as  computed  by  (12)  via  both  the  IMU  acceleration  and  angular  rate  values  -  are 
over-plotted  in  red.  These  clearly  show  the  decoupled  nature  of  the  motion,  where  there  is  good 
agreement  to  alternating  portions  of  the  acceleration  components.  The  three  other  separate 
contributions  to  the  acceleration  components  are  shown  over-plotted  in  Figure  22.  These 
reaffirm  the  notion  based  on  Figure  21  that  the  motion  is  decoupled. 

As  in  the  case  of  the  purely  rotational  motion  of  Data  Set  B,  it  is  instructive  to  compare  the 
velocity  components  in  the  initial  body  frame  as  derived  by  integrating  the  gravitationally- 
compensated  accelerometer  components  to  those  derived  using  the  IMU  angular  rates  via  (8). 
These  are  respectively  represented  as  the  black  and  red  curves  of  Figure  23.  As  before,  the  linear 
contribution  in  (8)  is  neglected,  not  because  it  is  non-existent,  however,  but  because  it  cannot  be 
determined  from  the  IMU  angular  rates. 

The  fact  that  the  linear  and  rotational  motions  are  uncoupled  in  Data  Set  C2  allows  us  to  readily 
correlate  the  black  and  red  curves  of  Figure  23.  This  provides  us  with  the  means  to  completely 
characterize,  and  therefore  correct  for,  the  drift  in  the  gravitationally-compensated  accelerometer 
components.  If  we  use  the  simple  approach  of  characterizing  the  mismatch  of  the  correlated  data 
by  a  polynomial,  we  find  that  the  choice  of  quadratics  are  sufficient  (refer  to  Figure  24).  It 
should  be  noted  that  the  mismatches  we  are  attempting  to  characterize  here  are,  via  (8),  the 
velocity  components  of  the  pivot  point  in  the  initial  body  frame  with  drifts  included. 
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Figure  17:  The  accelerometer-derived  (black)  vs.  the  gyroscope-derived  (green)  velocity  components  for  the 
IMU  in  the  initial  body  frame  for  Data  Set  B. 
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Figure  18:  The  accelerometer-derived  (black)  vs.  the  gyroscope-derived  (green)  IMU  positions  in  the  initial 
body  frame  for  Data  Set  B. 
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Figure  19:  The  gyroscope-derived  vs.  video  ground  truth  (red)  positions  for  sticker  4  in  the  video  frame. 
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Figure  20:  The  mismatches  between  the  gyroscope-derived  and  the  video  ground  truth  X  and  Y  positions  of 
Figure  19. 
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Figure  21:  Acceleration  components  as  measured  via  the  IMU  for  Data  Set  C2.  The  over-plotted  red  curves 
represent  the  linear  contribution  based  on  equation  (12). 
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Figure  22:  Acceleration  components  as  measured  via  the  IMU  for  Data  Set  C2.  The  over-plotted  colored 
curves  show  the  respective  contributions  to  the  acceleration  components. 
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Figure  23:  The  accelerometer-derived  (black)  vs.  the  gyroscope-derived  (red)  velocity  components  for  the 
IMU  in  the  initial  body  frame  for  Data  Set  C2.  The  gyroscope-derived  velocity  components  neglect  the  linear 
contributions  to  velocity,  as  a  matter  of  definition,  but  remain  much  more  stable. 
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The  quadratic-corrected  accelerometer-derived  velocity  components  are  represented  in  Figure  25 
by  the  green  curves.  It  is  obvious  that  although  the  agreement  with  the  red  curves  (i.e.  the 
velocity  components  derived  by  using  the  IMU  angular  rates)  has  improved,  the  match  remains 
far  from  perfect.  However,  judging  from  Figure  26,  the  positions  obtained  by  integrating  the 
quadratic-corrected  accelerometer-derived  velocity  components  (denoted  by  the  green  curves) 
show  reasonably  good  agreement  with  the  ground  truth  positions  extracted  from  the  video  data 
(denoted  by  the  red  curves),  especially  when  compared  to  the  positions  obtained  by  integrating 
the  uncompensated  accelerometer-derived  velocity  components  (denoted  by  the  black  curves).  In 
fact,  as  Figure  27  shows,  the  match  is  usually  within  a  few  centimeters  and,  as  Figure  28 
demonstrates,  adheres  generally  well  to  the  true  trajectory. 


Figure  24:  The  rms  fit  errors  for  polynomial  fits  of  given  order  to  the  mismatch  between  accelerometer- 
derived  and  the  gyroscope-derived  velocity  components  for  the  IMU  in  the  initial  body  frame  for  Data  Set 


C2. 


As  mentioned  previously,  the  means  to  completely  characterize,  and  therefore  correct  for,  the 
drift  in  the  gravitationally-compensated  accelerometer  components  is  at  our  disposal  here. 
However,  a  simple  polynomial  characterization  of  the  mismatch  appears  to  be  inadequate  for  the 
desired  level  of  positional  accuracy  that  we  know  is  achievable  (Figure  20).  The  question  then 
becomes  -  what  processing  steps  can  we  take  to  attain  the  highest  level  of  positional  accuracy? 
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Figure  25:  The  polynomial-corrected  accelerometer-derived  (green)  vs.  the  gyroscope-derived  (red) 
velocity  components  for  the  IMU  in  the  initial  body  frame  for  Data  Set  C2.  The  polynomials  chosen  for  the 
corrections  are  quadratics  for  all  three  components  (choice  based  on  Figure  24) 
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Figure  26:  The  video  frame  positions  obtained  by  integrating  the  uncompensated  accelerometer-derived 
velocity  components  (black)  vs.  those  obtained  by  integrating  the  quadratic-corrected  accelerometer-derived 
velocity  components  (green).  The  red  curves  represent  the  ground  truth  positions  extracted  from  the  video 
data  for  Data  Set  C2. 
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Figure  27:  The  mismatches  between  the  compensated  accelerometer-derived  and  the  video  ground  truth  X 
and  Y  positions  of  Figure  26. 
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3. 1.2.3.  Coupled  Rotational  and  Linear  Motion 

Data  Set  E  represents  the  case  where  the  rotational  and  linear  motions  of  the  body  are  coupled. 
The  x,  y  and  z  IMU  accelerometer  values  for  this  data  set  are  shown  in  Figure  29  and  Figure  30. 
In  Figure  29,  the  linear  contributions  -  as  computed  by  (12)  via  both  the  IMU  acceleration  and 
angular  rate  values  -  are  over-plotted  in  red.  These  clearly  show  the  coupled  nature  of  the 
motion,  where  there  is  now  generally  poor  agreement  between  the  curves.  The  three  other 
separate  contributions  to  the  acceleration  components  are  shown  over-plotted  in  Figure  30. 
These  reaffirm  the  notion  based  on  Figure  29  that  the  motion  is  coupled. 

As  before,  it  is  instructive  to  compare  the  velocity  components  in  the  initial  body  frame  as 
derived  by  integrating  the  gravitationally-compensated  accelerometer  components  to  those 
derived  using  the  IMU  angular  rates  via  (8).  These  are  respectively  represented  as  the  black  and 
red  curves  of  Figure  31.  Obviously,  the  fact  that  the  linear  and  rotational  motions  are  now 
coupled  makes  it  much  harder  to  correlate  the  two  curves.  However,  if  we  ignore  this  fact  and 
characterize  the  mismatch  for  data  based  on  the  previously  defined  criterion  (i.e.  the  intersection 
of  the  data  where  the  red  curves  of  vx  and  vy  exceed  a  threshold  and  are  within  the  time  interval 
of  the  collected  video  data)  by  quadratics  again,  we  obtain  the  quadratic-corrected 
accelerometer-derived  velocity  components  represented  in  Figure  32  by  the  green  curves. 


x  (m) 


Figure  28:  The  trajectory  of  sticker  4  based  on  the  ground  truth  positions  extracted  from  the  video  data 
(red)  vs.  the  trajectory  obtained  by  integrating  the  quadratic-corrected  accelerometer-derived  velocity 
components  (green)  for  Data  Set  C2. 
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Figure  29:  Acceleration  components  as  measured  via  the  IMU  for  Data  Set  E.  The  over-plotted  red  curves 
represent  the  linear  contribution  based  on  equation  (12).  Note  that  the  linear  contribution  is  not  decoupled 
from  the  rotational  contribution  (cf.  with  Figure  30). 
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Figure  30:  Acceleration  components  as  measured  via  the  IMU  for  Data  Set  E.  The  over-plotted  colored 
curves  show  the  respective  contributions  to  the  acceleration  components. 
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Figure  31:  The  accelerometer-derived  (black)  vs.  the  gyroscope-derived  (red)  velocity  components  for  the 
IMU  in  the  initial  body  frame  for  Data  Set  E.  The  gyroscope-derived  velocity  components  neglect  the  linear 
contributions  to  velocity,  as  a  matter  of  definition,  but  remain  much  more  stable. 
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Integrating  both  the  compensated  accelerometer-derived  velocity  components  (i.e.  the  green 
curves  in  Figure  32)  and  the  uncompensated  accelerometer-derived  velocity  components  (i.e.  the 
black  curves  in  Figure  31),  we  obtain  the  results  for  the  x  and  y  positions  in  the  video  frame 
shown  in  Figure  33.  The  red  curve  here  shows  the  ground  truth  positions  extracted  from  the 
video  data  for  Data  Set  E.  It  is  obvious  from  the  figure  that  there  is  a  significant  mismatch  in  the 
y  position  of  the  video  coordinate  system,  even  after  the  aforementioned  compensation.  What  is 
interesting  though  is  that  the  y  position  mismatch  is  to  first  order  linear  and  is  sufficiently 
characterized  by  a  third  order  polynomial  leading,  after  the  correction  process,  to  cm-scale 
agreements  in  positioning  demonstrated  by  Figure  34  and  Figure  35. 
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Figure  32:  The  quadratic-corrected  accelerometer-derived  (green)  vs.  the  gyroscope-derived  (red)  velocity 
components  for  the  IMU  in  the  initial  body  frame  for  Data  Set  E.  Note  that  only  the  IMU  data  in  the  time 
interval  of  Data  Set  E  (i.e.  0  -  17  s)  was  chosen  for  the  correction  process. 
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Figure  33:  The  video  frame  positions  obtained  by  integrating  the  uncompensated  accelerometer-derived 
velocity  components  (black)  vs.  those  obtained  by  integrating  the  quadratic-corrected  accelerometer-derived 
velocity  components  (green).  The  red  curves  represent  the  ground  truth  positions  extracted  from  the  video 
data  for  Data  Set  E. 
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Figure  34:  The  trajectory  of  sticker  4  based  on  the  ground  truth  positions  extracted  from  the  video  data 
(red)  vs.  the  trajectory  obtained  by  fully  compensating  as  described  in  the  text  (green)  for  Data  Set  E. 
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Figure  35:  The  mismatches  between  the  red  and  green  X  and  Y  positions  of  Figure  34. 
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3.2.  Electromagnetic  Induction  Sensor 

As  stated  in  the  introduction,  the  ultimate  goal  of  this  project  is  to  develop  a  handheld  sensor  that 
characterizes  a  target  for  discrimination  as  either  UXO  or  clutter,  and  that  does  this  with  a  low 
probability  of  false  alann.  In  order  to  accomplish  this,  a  crucial  ingredient  is  a  reliable  inversion 
algorithm  that  takes  the  electromagnetic  induction  (EMI)  data  and  backs  out  a  minimum  number 
of  parameters  to  make  the  decision.  These  target-characterizing  parameters  model  the  very 
specific  spatial  pattern  of  the  EMI  data  the  target  produces  ([6]— [7]).  As  a  consequence,  the 
inversion  -  which  incorporates  this  model  -  requires  the  data  to  be  spatially  well-located  for 
proper  characterization.  In  fact,  it  was  determined  in  ESTCP  Project  200108  (EM61-HH  with 
template)  that  position  accuracies  of  a  couple  of  cm  were  required  to  perfonn  dependable 
inversions.  This  was  based  on  the  attainment  of  squared  correlations  (between  the  data  and  the 
model)  greater  than  0.99  for  compact  targets. 

3.2.1.  Inversion  Results 

If  we  take  the  EMI  and  compensated  IMU-derived  trajectory  of  Figure  33  for  Data  Set  E  and 
carry  out  the  inversion  (refer  to  [6]  for  specifics),  we  get  the  result  shown  in  Figure  36.  The 
resulting  target  parameters  are  comparable  to  those  obtained  using  the  ground  truth  positions 
extracted  from  the  video  data.  Both  inversion  results  are  presented  in  Table  3  and  Table  4  below. 
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(°) 

Rz 
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-0.076 

-0.107 

-0.367 

0.2 

-0.3 

111.4 

0.989 

Video 

-0.076 

-0.107 

-0.354 

0.0 

-0.5 

114.7 

0.991 

Table  3:  The  target  position,  orientation  and  squared  correlation  inversion  results,  for  the  two  cases  where 
the  IMU-derived  trajectory  and  the  video  ground  truth  trajectory  were  used  to  position  the  EMI  data  of 
Data  Set  E. 
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0.44 

0.18 

0.48 

0.23 

0.09 

0.32 

0.12 

0.03 

0.19 

0.06 

0.02 

Video 

0.71 

0.38 

0.19 

0.40 

0.19 

0.09 

0.27 

0.10 

0.04 

0.16 

0.05 

0.02 

Table  4:  The  responses  (units  of  m3)  that  are  induced  when  the  primary  field  is  aligned  in  turn  along  each 
of  the  three  principal  axes  of  the  target,  for  each  time  gate  of  the  EM61-HH. 


From  Table  3,  we  see  that  the  squared  correlations  are  on  the  lower  end  of  the  range  required  for 
reliable  characterization.  Although  the  position  and  orientation  results  for  the  target  remained 
relatively  consistent,  both  fits  predicted  the  target  to  be  a  few  cm  deeper  than  in  reality.  The 
target  here  is  a  horizontally  emplaced  40  mm  projectile  32  cm  below  the  sensor  head  (refer  to 
Figure  3).  In  addition,  the  resulting  target  response  coefficients  given  in  Table  4  (which  are 
strongly  correlated  with  the  dimensions  of  the  target,  and  so  should  yield  Py  ~  (3Z)  demonstrated 
much  more  spread  in  the  transverse  components  than  would  normally  be  desired  for  reliable 
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Figure  36:  Contours  of  the  EM61-HH  data  (blue  curves)  vs.  those  of  the  model  based  on  the  inversion 
target  parameter  results  (red  curves)  for  Data  Set  E  with  the  compensated  IMU-derived  positions  (black 
dots)  employed.  The  grey  dots  represent  the  ground  truth  positions  extracted  from  the  video  data  (for 
reference). 
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characterization.  Similar  disappointing  results  were  observed  when  inverting  the  EMI  data  of 
Data  Set  C. 


Significantly,  both  the  video  ground-truthed  and  IMU-positioned  inversions  have  fit  coherence 
substantially  poorer  that  typically  found  with  the  template  data  obtained  in  ESTCP  Project 
200108.  The  fit  coherence  (R2)  is  the  square  of  the  correlation  coefficient  between  the  data  and 
the  dipole  model  fit  to  the  data.  a/(1-R2)  is  a  measure  of  the  mismatch  between  the  model  fit  and 
the  data.  Template  data  is  typically  found  to  have  fit  coherences  of  about  0.999.  A  significant 
difference  between  these  data  and  the  template  data  is  that  these  data  were  obtained  in  a  dynamic 
environment  where  the  sensor  was  in  motion  rather  than  sitting  still  on  a  grid  point.  This 
suggests  the  need  for  a  closer  look  at  the  EM61-HH  response  in  a  dynamic  mode  of  operation, 
and  incorporating  this  into  our  current  model.  This  is  treated  in  the  following  section. 
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4.  EM61-HH  Dynamic  Response 


When  sweeping  the  EM61-HH  back  and  forth  over  a  target,  it  is  noted  that  the  peak  response  is 
both  delayed  and  distorted  relative  to  the  actual  target  location.  The  receiver  output  of  this  sensor 
is  analog  integrated  with  a  filter  that  both  shifts  and  distorts  the  sensor’s  response.  This  filter  can 
sufficiently  distort  signal  shape  to  limit  the  inversion  process.  This  has  a  significant  impact  on 
the  accuracy  of  target  parameter  estimates  from  dipole  inversion  at  sweeping  speeds  of  more 
than  a  few  cm/sec.  At  such  speeds,  it  is  necessary  to  account  for  the  sensor's  temporal  response. 

4.1.  Dynamic  Response  Characteristics 

The  manufacturer  (Geonics,  Ltd.)  informed  us  in  private  correspondence  that  EM61-HH  system 
can  be  represented  as  an  underdamped  second-order  system  given  by  the  standard  transfer 
function  of  the  form 


(13) 


H(s)  = 


co„ 


s 2  +  2  Ay  s  +  co  ~ 

~  n  n 


where  con  is  the  natural  frequency,  with  con>  0 ;  and  C  is  the  damping  factor,  with  0  <  C,  <  1 . 
To  determine  the  applicable  values  for  con  and  (T  for  our  sensor,  we  turned  to  measurements. 
Using  a  small  wire  coil  that  switches  between  high  and  low  impedance  as  triggered  by  a  pulse 
generator,  the  EM61-HH  response  to  a  five  second  on  /  five  second  off  signal  was  measured 
(Figure  37).  The  analytical  expression  for  the  step  response  of  such  a  system  is 


(14) 


y( 0  =  1  -  r - —  sin(<y  Vw2 1  +  cos-1 O 


and  was  found  to  fit  the  measured  data  well  with  con  =  6.2  and  C  =  0-78  (Figure  38). 
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Figure  37:  The  measured  EM61-HH  response  to  a  sequence  of  five  second  on  /  five  second  off  pulses. 


Figure  38:  The  normalized  EM61-HH  response  to  a  five  second  on  /  five  second  off  signal.  Data  from  all  the 
pulses  in  Figure  20  are  used.  The  green  curve  represents  the  best  fit  to  the  data  using  (13)  with  coefficients 

a>„  =6.2  and  ^  =  0.78. 
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4.2.  EM61-HH  Two-Dimensional  Sweep  Data 


The  original  apparatus  used  in  the  2D  data  collection  process  was  modified  in  several  ways  to 
minimize  any  extraneous  sources  of  error. 

First,  the  equipment  (i.e.  the  Crossbow  IMU400  and  the  Geonics  EM61-HFI  system  -  electronics 
&  data  logger  included)  was  attached  to  a  more  rigid  flat  board  and  allowed  to  slide  on  a  newly 
designed,  level  and  stable,  platform  (Figure  39).  The  EM61-HFI  coil  head  was  detached  from  its 
pole  and  placed  at  one  end  of  the  board,  while  the  EM61-FIFI  electronics  &  data  logger  along 
with  the  IMU,  were  placed  at  the  other  end  of  the  board  (Figure  40).  As  before,  movement  of  the 
board  was  constrained  to  rotate  about  a  pivot  point,  and  the  pivot  point  constrained  to  translate 
along  a  linear  groove  built  into  part  of  the  platform,  allowing  for  a  variety  of  simple  planar 
motions  containing  both  translational  and  rotational  contributions.  This  time,  however,  the  IMU 
was  centered  on  the  pivot  point.  Again,  the  video  camera  was  located  directly  above  the 
platform,  pointing  downwards. 


Figure  39:  New  overall  setup  for  2D  motion  measurements 


Second,  three  bright  LED’s  were  placed  on  the  EM61-HH  coil  head  to  enable  the  automated 
extraction  of  the  ground  truth  positioning  information  from  the  video  data  (Figures  40  and  41). 
To  aid  in  the  latter,  a  tarp  was  used  to  shade  the  platform  and  the  video  camera  exposure  was 
adjusted  to  reduce  the  background  (i.e.  non-LED)  intensities. 
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Figure  40:-  Close  up  of  rigid  flat  board 


Finally,  a  peg  board  with  holes  spaced  1”  apart  in  a  regular  2D  grid  pattern  (Figure  41)  was 
imaged  by  the  camera  at  the  height  of  the  LED ’s.  This  was  used  to  find  a  polynomial  fit  that 
would  map  camera  pixel  locations  to  a  calibrated  rectangular  X-Y  coordinate  system.  The 
location  of  each  LED  was  calculated  from  the  centroid  of  its  image  over  an  eight  pixel  square. 
When  the  board  was  moved  quickly  (1  m/s),  the  image  of  each  LED  would  noticeably  smear  out. 
As  a  sanity  check,  the  separation  of  the  two  outer  LED’s  was  measured  at  25.4  cm  and  compared 
to  the  mapped  images.  The  mean  separation  noted  in  the  mapped  images  was  25.3  cm  with  a 
standard  deviation  of  0.2  cm  as  the  board  was  swept  about. 

The  IMU  data  was  initially  collected  on  an  IPAQ  PDA  (as  shown  in  Figure  40),  but  the  program 
to  do  this  turned  out  to  be  unstable.  Instead,  a  laptop  with  the  GyroView  program  provided  by 
the  manufacturer  was  used.  This  meant  that  a  serial  cable  (from  IMU  to  laptop)  had  to  be 
carefully  dragged  around  while  the  board  was  being  moved. 

The  EMI  and  IMU  data  were  acquired  using  independent  clocks  and  then  manually  synchronized 
after  the  fact.  In  principle,  this  should  introduce  one  extra  unknown  representing  the  relative  time 
lag  between  the  two  data  streams.  However,  further  complications  arose  due  to  limitations 
associated  with  the  EM61-HH  data  logger.  Specifically,  sampling  rate  stability  issues  combined 
with  coarse  resolution  of  the  data  logger  clock  conspired  to  blur  the  true  relative  time 
relationship  between  the  two  data  streams.  In  order  to  address  this  problem,  by  the  summer  of 
2005,  both  EMI  and  IMU  data  were  being  collected  with  a  common  time  stamp  on  the  laptop 
using  a  LabView  DAQ  routine. 
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Figure  41:  Automatic  extraction  of  positioning  ground  truth  from  video  frame  data  via  intensity  threshold 
using  LEDs 


4.3.  Inversion  Results 

The  transfer  function  with  the  established  values  for  our  sensor  was  incorporated  into  the  dipole 
inversion  model.  Figure  42  shows  comparisons  of  the  measured  EM61-HFI  signal  during  a 
sweep  over  a  target  and  dipole  fits  using  the  static  and  the  dynamic  response  models.  The  data 
are  from  three  passes  near  the  middle  of  the  sensor  trajectory.  The  dynamic  response  model 
clearly  does  a  better  job  of  representing  the  observed  response. 
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Figure  42:  Comparison  of  static  and  dynamic  response  models  with  measured  response. 
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Figure  43  compares  inversion  results  with  the  static  (upper  plot)  and  dynamic  (lower  plot) 
response  models.  The  data  were  collected  on  a  test  stand  with  the  EM61-HH  swept  back  and 
forth  over  the  target.  In  these  figures,  the  primary  (largest)  P  is  plotted  along  the  horizontal  axis, 
and  the  secondary  Ps  are  plotted  along  the  vertical  axis,  with  a  symbol  at  the  mean  value  and  a 
line  running  between  the  two.  For  the  test  objects  (two  spheres  and  two  projectiles),  the  two 
transverse  beta  values  are  expected  to  be  equal.  Furthermore,  for  the  two  spheres  considered,  all 
three  betas  (i.e.  longitudinal  and  transverse)  are  expected  to  be  equal.  The  dotted  line  represents 
the  line  where  longitudinal  and  transverse  values  are  equal.  The  dynamic  model  generally 
produces  tighter  eigenvalue  clusters  and  better  shape  characterization  (secondary  Ps  more  nearly 
equal)  than  does  the  static  model. 
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Figure  43:  Inversion  results  showing  ps  when  the  dynamic  response  is  not  taken  into  account  (upper 
panel)  versus  when  it  is  taken  into  account  (lower  panel). 
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5.  Arcsecond  Laser-Based  Positioning 


Position  ground  truth  becomes  a  significant  issue  when  we  consider  three-dimensional  motion. 
The  video  system  is  only  appropriate  for  2D  motions.  A  laser  positioning  system  developed  by 
Arcsecond  (subsequently  acquired  by  Metris)  was  adapted  for  UXO  field  surveys  in  ESTCP 
project  MM-0129.  Arcsecond  refers  to  their  systems  as  “Indoor-GPS”  positioning  systems. 
Much  of  what  they  do  is  used  for  very  precise  measurements  in  manufacturing  settings.  The 
system  is  based  on  fixed  rotating  laser  transmitters  that  they  liken  to  GPS  satellites.  The 
transmitters  have  two  rotating  laser  fans  that  are  tilted  relative  to  each  other,  plus  a  timing  light 
strobe.  By  measuring  the  timing  of  the  two  laser  fans  relative  to  the  timing  strobe,  one’s  relative 
angular  position  from  a  transmitter  can  be  measured.  Given  a  set  of  two  or  more  transmitters  at 
known  locations,  all  of  these  angular  measurements  can  be  used  to  calculate  one’s  position.  The 
ESTCP  MM-0129  system  makes  use  of  four  rigidly  attached  light  sensors  in  a  tetrahedral 
configuration  to  measure  both  the  survey  platfonn’s  position  and  orientation.  In  post-processing 
software,  the  displacements  from  the  Arcsecond  sensors  to  the  EMI  sensor  can  be  used  to 
calculate  the  3D  position  and  orientation  of  the  EMI  coil.  Limited  by  the  laser  rotation  rate 
around  40Hz,  the  position  of  the  platform  is  measured  at  a  rate  of  20  samples  per  second. 

5.1.  Test  Setup 

Our  testing  of  the  Arcsecond  system  for  this  project  was  reported  in  [8].  The  tests  were  done  on  a 
level  wooden  platform.  Test  objects  could  be  placed  roughly  30  centimeters  below  the  platform. 
Photos  of  the  test  setups  are  shown  in  Figure  44.  All  of  the  data  presented  here  is  over  a  four 
inch  diameter  carbon  steel  sphere.  The  2D  video  and  2D  Arcsecond  position  data  were  collected 
simultaneously.  EM61-HH  data  was  collected  directly  onto  a  laptop  PC  at  a  rate  of  15  samples 
per  second.  The  PC  was  also  used  to  collect  the  Crossbow  IMU  data. 

The  2D  tests  were  done  by  attaching  the  equipment  to  a  flat  board  and  sliding  the  board  along 
the  platform  surface.  The  EM61-HH  coil  head  was  detached  from  the  usual  pole  and  located  at 
one  end  of  the  board.  At  the  other  end  of  the  board,  the  Geonics  electronics,  the  Crossbow  IMU, 
and  the  Arcsecond  sensors  were  attached  (Figure  44a  and  44b).  As  noted  previously,  the  sensor 
location  was  tracked  using  video  imaging  of  LEDs  on  the  coil.  The  camera  frame  rate  was  30 
frames  per  second. 

The  Arcsecond  system  used  for  the  2D  tests  had  a  tetrahedral  configuration  of  four  sensors 
attached  to  the  sensor  board  75  cm  back  from  the  coil.  The  displacements  from  each  Arcsecond 
sensor  to  the  EM  coil  were  measured  for  post-processing  the  coil’s  position.  When  stationary, 
the  RMS  noise  in  the  sensor  positions  is  on  the  order  of  0.2  millimeters  and  the  calculated  EM 
coil  position  noise  is  around  0.6  mm.  For  the  2D  tests,  only  two  laser  transmitters  were  used. 
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Figure  44:  2D  and  3D  EM61-HH  sensor  positioning  experiments.  In  (a)  sideview  and  (b)  video 
camera  view,  the  2D  sensor  board  is  shown  with  EM61  coil,  LEDs,  Arcsecond  sensors  and  Crossbow 
IMU.  In  (c)  the  3D  setup  is  shown  with  standard  EM61-HH  pole,  Arcsecond  sensors  and  IMU. 


55 


For  the  3D  test,  the  EM61-HH  coil  was  reattached  to  the  standard  Geonics  rig.  A  flat  board  and  a 
surveyor  pole  were  attached  to  the  rig  to  hold  both  the  Crossbow  IMU  and  the  Arcsecond 
sensors.  A  shoulder  strap  held  the  rig  over  the  experimenter’s  shoulder,  and  the  EM  coil  could 
be  swept  side-to-side  in  a  typical  field  survey  interrogation  of  an  object  (see  Figure  44c).  For 
each  data  set,  the  coil  would  start  flat  and  be  lifted  from  the  test  platform  for  the  sweep  and 
placed  flat  on  the  platform  surface  at  the  end  of  each  sweep. 

5.2.  Test  Results 

Figure  45  plots  sample  trajectories  of  the  EM  coil  from  the  three  positioning  data  sets.  An 
example  of  the  concurrently  collected  video  and  Arcsecond  positions  are  shown  in  45a  and  45c. 
The  video  data  (red  curve)  was  shifted  in  X,  Y,  and  time  to  overlay  with  the  Arcsecond 
positions.  To  compare  the  two  sets,  the  video  data  was  interpolated  to  the  Arcsecond  times. 
Across  the  central  region  of  the  video  data,  the  standard  deviations  in  the  differences  in  X  and  Y 
were  between  2  and  4  mm.  The  mean  horizontal  distance  between  the  two  position 
measurements  was  3.5  mm  with  a  standard  deviation  of  2.3  mm.  Larger  differences  were  noted 
around  the  edges  of  the  video  image.  There  is  a  distinct  trend  to  differences  at  the  edge  and  there 
maybe  a  residual  error  to  the  applied  lens  correction.  The  video  data  assumes  a  flat  board,  but  the 
Arcsecond  data  in  45c  shows  almost  2  cm  in  variation.  Along  the  Y  direction,  there  is  a  distinct 
slope  to  the  measured  Z  on  the  order  of  1.0  cm  over  a  distance  of  100  cm.  Over  a  narrow  central 
region  of  the  sweeps,  the  standard  deviation  in  Z  is  2.6  mm.  Based  on  the  reasonable  agreement 
between  these  two  systems  we  conclude  that  their  positioning  is  probably  good  to  5  mm  or 
better.  Figure  45b  and  45d  plot  a  sample  3D  trajectory  of  the  EM  coil.  The  measured  variation  in 
Z  as  the  coil  is  lifted  and  swept  side-to-side  over  the  platform  is  on  the  order  of  10  cm. 

The  EMI  data  were  fit  with  the  standard  induced  magnetic  dipole  response  model,  oincluding  the 
dynamic  response  of  the  EM61-HFI.  The  model  parameters  are  object  position  (x,  y,  z),  object 
orientation  (6,  (f>,  <//),  and  the  object’s  magnetic  polarization  responses  along  its  primary  axes 
(/?/,  J32,  J3j).  Because  the  test  object  was  a  sphere,  we  simplified  our  fits  with  a  single  J3  value. 
Another  important  aspect  of  the  inversion  model  was  accounting  for  the  temporal  response  of  the 
EM61-HH 

Results  of  fitting  the  measured  data  with  the  different  positioning  systems  are  shown  in  Figures 
46  and  47.  The  black  symbol/curves  are  the  measured  data  and  the  green  curve  is  the  best  model 
fit.  In  Figure  46,  the  data  is  plotted  as  a  function  of  time.  In  Figure  47,  it  is  plotted  as  a  function 
of  the  X  position.  Note  the  displacements  in  X  of  the  peak  signals  due  to  the  sensor  response;  the 
actual  sphere  location  is  between  the  peaks. 

For  the  2D  video/ Arcsecond  data,  the  measurements  over  the  sphere  were  repeated  and  inverted 
seven  times.  For  the  3D  data,  seventeen  sets  of  measurements  were  made.  Figure  48  presents  the 
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Figure  45:  Sample  X-Y  and  X-Z  plots  of  the  EM61-HH  coil  trajectories.  In  (a)  and  (c),  the  2D  coil 
motion  on  the  sensor  board  as  measured  by  the  Arcsecond  (black  symbols)  and  the  video  (red  curve). 
In  (b)  and  (d),  the  3D  motion  of  the  coil  on  the  standard  EM61-HH  pole.  The  green  symbol  is  the 
object  location. 


inversion  results  for  all  of  this  data  (red  X  -  video  fit,  green  triangle  -  2D  Arcsecond,  black 
diamond  -  3D  Arcsecond). 

The  fitted  depths  and  response  P s  are  plotted  in  the  top  panel  of  Figure  48.  These  two  fit 
parameters  are  highly  correlated  in  the  inversion  model.  Given  measurements  with  both  signal 
noise  and  positioning  error,  the  inversion  model  tends  to  vary  these  two  parameters  the  most  in 
an  attempt  to  get  the  best  model  fit  to  the  data.  For  the  three  positioning  sets,  the  mean  sphere 
P s  are:  6.70  -  2D  Arcsecond,  6.88  -  2D  video,  and  6.94  -  3D  Arcsecond.  The  standard 
deviations  are:  0.18,  0.14,  and  0.27,  respectively.  The  sphere’s  depth  was  not  noted  accurately 
enough  to  check  on  the  fitted 
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Figure  46:  Time  rasters  of  EM61-HH  data  (black  curve  and  symbols)  fitted  to  model  (green  curve). 
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Figure  47:  Plots  along  X  of  EM61-HH  data  (black  curve  and  symbols)  fitted  to  model  (green  curve). 
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Figure  48:  Results  of  fitting  EM61-HH  data  with  model  using  different  positioning  data. 


depths.  It  was  roughly  29  cm  below  the  center  of  the  EM  coil  when  it  was  resting  flat  on  the 
platform.  For  the  three  coordinate  systems,  z  =  0  was  only  roughly  matched  to  this  coil  height. 
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In  the  lower  panel  of  Figure  48,  a  measure  of  the  fit  quality  is  plotted  versus  the  average  swing 
speed  of  the  EM61-HH  coil.  The  peak  speeds  were  on  the  order  of  2  to  3  times  greater.  The  fit 
quality  parameter  is  the  square  of  the  coherence  between  the  best  model  fit  and  the  data.  There  is 
a  clear  trend  to  poorer  fits  as  the  coil  velocity  increases.  This  is  not  surprising  given  the 
positioning  system  limitations  of  video  frame  rate  (30  fps)  and  laser  rotation  rate  (40Hz).  At  a 
speed  of  1  m/s,  the  coil  has  moved  3.3  cm  in  a  video  frame  and  2.5  cm  during  a  single  laser  spin. 
Man  portable  survey  systems  move  at  comparable  speeds  and  vehicle  survey  platforms  exceed  it. 

In  an  effort  to  simulate  the  error  in  inverting  the  polarization  responses  as  a  function  of 
positioning  errors,  a  simple  Monte  Carlo  simulation  was  run.  Static  measurements  on  a  fixed 
grid  with  exact  positioning  were  taken  over  a  40mm  projectile.  The  data  were  inverted  for  the  // 
responses.  The  same  data  was  inverted  repeatedly  with  increasing  random  errors  added  to  the 
known  grid  positions.  The  results  are  displayed  in  Figure  49,  which  shows  the  variability  in  the 
fitted  polarization  responses  versus  RMS  position  error  in  the  top  panel  and  the  coherence 
squared  of  the  fit  the  lower  panel.  At  the  level  of  3  mm  RMS  position  error,  the  average  fit 
quality  is  about  0.998  and  the  spread  in  primary  response  is  15%.  This  is  consistent  with  the 
Arcsecond  3D  fits  having  an  average  fit  quality  of  0.9982  and  a  P  spread  of  13%  (-0.9/6.94). 

To  summarize,  in  the  2D  tests,  the  positions  mapped  were  found  to  match  video  mapped  images 
at  RMS  position  differences  of  3-4  mm.  Inverted  2D  and  3D  positioned  data  had  high  fit 
qualities  of  0.9982.  The  fitted  magnetic  polarization  responses  had  a  narrow  range  of  values 
(13%)  over  a  set  of  seventeen  3D  measurements  of  a  steel  sphere.  However,  the  fit  quality 
decreased  at  faster  sensor  motion  rates,  indicating  larger  positioning  errors  as  speeds  approach  1 
m/s.  Simulations  of  EMI  data  inversion  with  randomized  position  errors  came  up  with  similar  fit 
qualities  and  response  parameter  errors  at  an  RMS  position  error  level  of  3mm.  All  indications 
are  that  this  laser  based  system  is  tracking  the  EMI  sensor  with  sub-centimeter  accuracy,  and  is 
suitable  for  ground-truthing  3D  IMU  position  data. 
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Figure  49:  Monte  Carlo  simulation  of  inverting  40mm  projectile  data  with  randomly  generated 
position  errors.  Dashed  line  at  RMS  position  error  of  3  mm  with  an  average  lit  coherence  of  0.998 
and  a  primary  b  spread  of  15%  (0.05/0.34). 
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6.  Three  Dimensional  Motion 


With  general  three-dimensional  motions  it  is  necessary  to  have  some  sort  of  external  aid  in  order 
to  properly  compensate  the  IMU  measured  quantities  and  bound  the  random  walk  errors.  This 
may  involve  either  deploying  additional  sensors  to  get  periodic  independent  confirmation  of  the 
true  or  derived  quantities,  or  periodically  forcing  the  IMU  to  be  in  situations  where  the  true  or 
derived  quantities  are  known.  The  simplest  scheme  is  one  where  zero  velocity  updates  are 
periodically  taken. 

6.1.  Zero  Velocity  Update  Method 

The  quantities  measured  by  the  IMU  are  the  angular  rate  vector,  65 ,  and  the  acceleration  vector, 
a  .  The  angular  rate  can  be  expressed  in  the  IMU  body  frame  of  reference  as 

(15)  S(0  =  s„ -©(O  +  MO  +  nJO 


where  or(/)  represents  the  true  angular  rate  vector;  b(i  (/)  and  n(i  (/)  represent,  respectively,  the 
bias  and  zero-mean  noise  vectors  in  the  angular  rate  measurements.  A  constant  calibration  scale 
factor  sm  is  included  for  generality.  Comparisons  with  ground  truth  showed  that  the  calibration 
factor  can  be  set  to  1  for  our  application.  The  acceleration  vector  is 

(16)  a(0  =  s,  •a(i)  +  Rc^(i)g  +  ba(i)  +  na(0 


where  a (t)  represents  the  true  acceleration  vector,  free  of  the  influence  of  gravity.  ba  (/)  and 
na(t)  represent,  respectively,  the  bias  and  zero-mean  noise  vectors  in  the  acceleration 
measurements;  and  sa  represents  a  calibration  scale  factor  vector,  equal  to  1  for  our  application. 
Re_>i(0  is  the  rotation  matrix  that  transforms  earth’s  gravitational  acceleration  vector  g  from  the 
earth-fixed  frame  to  the  IMU  body  frame  at  time  t. 


(17)  R,_J0 


cos//(t)cos9(/j 

-  s  in yAt)  c  o  s$f ) + c  o  s jAt)  s  in 6{t)  sin (f{t) 
sin^i)sin$Y) + cos^f)sin$(f)cos$/) 


sin(//(/)cos<9(t) 

c  o  s/A,t)  c  o  s$f) + s  im/Aj)  s  in$(t)  s  mjit) 
-  c  o  ^/At)  s  in$Y ) + s  in yAt)  s  m6(i)  c  o  s$7) 


-sirkf/j 

cos$(f)sin$Y) 

cos$(t)cos$/) 


where  y/(t ) ,  0(t)  and  rj>(t)  are  the  yaw,  pitch  and  roll  angles,  respectively  (see  Appendix  A). 
Hence 


"  0" 

-  sin  0{t) 

(18) 

Re^(Og  =  Re_>*(0 

0 

= 

cos  0{t)  sin  <j)(t ) 

cos  0{t)  COS  (j){t) 

The  trajectory  of  a  moving  IMU  body  in  the  local  earth-fixed  frame  is  determined  by  isolating 
the  true  acceleration  vector  a(/)  of  (16)  and  transforming  this  to  the  local  earth- fixed  frame, 
before  integrating  the  quantity  twice  -  once  to  obtain  the  velocity  vector,  and  a  second  time  to 
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obtain  the  position  vector.  The  task  of  isolating  a(t)  from  the  measured  acceleration  vector, 
however,  demands  an  accurate  instantaneous  accounting  of  both  the  gravitational  contribution  to 
acceleration  and  the  bias  vector,  as  well  as  knowledge  of  the  scale  factor  vector. 

In  order  to  properly  account  for  the  gravitational  contribution  to  acceleration,  the  evolving 
orientation  of  the  IMU  body  frame  relative  to  the  earth-fixed  frame  is  necessary.  The  orientation 
time  series  is  computed  by  integrating  the  true  angular  rate  vector  io(7)  of  (15)  to  obtain  the 
Euler  yaw,  pitch  and  roll  angles.  Repeating  (4),  this  represents  solving  a  set  of  coupled 
differential  equations 


(19) 


l 

sin^(t)  tan<9(t) 

cos  </>(t)  tan  9(t) 

0{t) 

= 

0 

cos  <f>(t) 

-  sin  <f>(t ) 

y{t)_ 

0 

sin  <f>(t) 

cos^(t) 

cos  6{t) 

cos6*(t) 

®,(0 

®y(i) 

®-(0 


where  y/(t) ,  0{t)  and  are  the  yaw,  pitch  and  roll  angles,  respectively. 

As  with  the  acceleration  vector,  the  task  of  isolating  (o(t)  from  the  measured  angular  rate  vector 
also  demands  an  accurate  instantaneous  accounting  of  the  bias  vector,  as  well  as  knowledge  of 
the  scale  factor  vector. 

It  should  be  pointed  out  that  even  if  one  were  able  to  perfectly  compensate  for  s(0 ,  b(i|  (/) ,  sa , 
R_*(0g  and  ba(t)  from  the  measured  vector  IMU  quantities,  boundlessly  growing  errors 
would  inevitably  still  arise  due  to  the  integration  operations.  These  errors,  known  as  random 
walk  errors,  occur  as  a  result  of  integrating  noise  terms,  whereby  the  variance  of  the  estimates 
increases  proportionately  to  the  integration  time. 

During  a  zero  velocity  update,  c o(t)  =  0  and  a (t)  =  0,  so  that  we  have,  respectively,  from  (15)  and 
(16), 

(20)  fi(t)  =  b(0(t)  +  n(0(t) 

and 

(2 1 )  a(0  =  RZb  g  +  b,  (0  +  n,  (0 

The  term  R  Zb  g  in  (21)  denotes  an  orientation-dependent  bias  vector  that  is  constant  for  the 
duration  of  the  zero  velocity  update,  and  is  represented  explicitly  via  (18)  with  9{t)  and  </>(t)  as 
constants.  Since  to(t)  and  a (t)  are,  to  first  order,  well  represented  by  linear  regression  vectors 
during  a  zero  velocity  update,  we  choose  to  write 

(22)  b„(0  =  Aa  +BJ  and  ba(f)  =  Aa+Baf 

Making  use  of  (22),  the  true  IMU  quantities  of  (15)  and  (16)  can  be  expressed,  in  the  mean,  as 
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(23) 


® (0  =  (l  +  8»)-  (®(0  -  (A«  +  B^)) 


and 

(24)  a«  =  (1  +  8,).  (a«)  -  R„„«g  -  (A,  +  B.0) 

where  1  +  5 (0  =sm_1;  l  +  5a  =sa_1;  1  is  the  identity  vector;  and  50) ,  5a  are  scale  factor  error 
vectors.  Furthermore,  transforming  (24)  from  the  IMU  body  frame  to  the  earth-fixed  frame  gives 

(25)  a *(<)  s  R  (/)a(i)=  R^(<)[(l  +  «.)■  (W)  ~  (A,  +  B.i))]-  (l  +  8,)  g 

where  R4^e(0  =  Rr  ;  (t)  is  the  transpose  of  (17),  and  Rr  (t)Re^.j(t)  =  I ,  the  identity  matrix. 

From  (23),  it  is  clear  that  nine  parameters  are  necessary  to  fully  specify  to(t) .  An  additional  six 
parameters  are  needed  to  accurately  determine  the  yaw,  pitch  and  roll  angles  from  t o(t)  via  (19); 
three  related  to  the  initial  angles,  and  three  related  to  the  final  angles.  The  former  are  crucial  for 
proper  initialization  of  (19),  while  the  latter  are  required  to  confine  the  random  walk  errors  that 
arise  due  to  the  integration  operations.  An  appropriate  choice  for  these  parameters  are  the  error 

angles  8y/i  =  y/i  -  y/ 1 ,  80,  =  6j  -  0. ,  8<l>i  =  -  fa ,  S\//  f  =  yz  f  -  y/  f ,  80 f  =  0  f  -  0  f  and 

8(/)f  =  (f)  f  -  <p  f ,  where  y/l ,  6*,  and  fa  are  the  true  initial  yaw,  pitch  and  roll  angles,  respectively; 
i/7 ; ,  0 1  and  fa  are  the  estimated  initial  yaw,  pitch  and  roll  angles,  respectively;  and  y/  f ,  0  f , 
<j)f ,  and  fa  f ,  0 f  ,  <j>f  are  the  true  and  estimated  final  angles,  respectively. 

By  starting  and  ending  with  zero  velocity  updates,  both  the  initial  and  final  pitch  and  roll  angles 
can  be  estimated  from  the  measured  acceleration  vector  using  (18),  (21)  and  (22).  That  is, 

-  sin  0ZV 

R«Zlg-  COS 9”  sin/-  g  =  »(0-(A.+B.O 

COS0ZI  COS0ZI 

which  leads  to 


(26)  0  =  sin-  - 


.  J  (a x( t)-(Aax+BaJ)) 


and  </)  =  sin- 


.  J(“yW~(Aay  +Baj)) 


g  COS0' 


If  it  is  assumed  that  estimates  for  the  initial  and  final  yaw  angles  are  also  available,  then,  by 
specifying  the  six  angle  error  parameters,  the  true  initial  and  final  angles  are  fully  determined. 
This  allows  needed  linear  corrections  to  be  made  to  each  of  the  y/(t) ,  0{t)  and  <j>(t)  solutions  to 
(19). 

From  (25),  to  fully  detennine  a e(t)  then,  we  need  the  15  parameters  just  described  to  completely 
specify  the  rotation  matrix,  and  an  additional  9  parameters  for  a  total  of  24  parameters.  The 
velocity  vector  -  which  is  the  integration  of  ae(t)  -  is  consequently  fully  determined  since  we 
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know  what  the  true  initial  and  final  velocities  are,  thereby  allowing  linear  corrections  to  be 
applied.  If  we  start  and  end  at  the  same  position,  the  same  can  be  said  of  the  position  vector  - 
which  is  the  integration  of  the  velocity  vector. 

With  all  of  the  assumptions  made  above,  we  can  therefore  completely  specify  the  position  vector 
in  the  earth- fixed  frame  between  two  co-located  zero  velocity  updates  with  24  parameters: 

K’  56 n  80  f  and  S(j)  f 

6.2.  Test  Setup 

A  considerable  simplification  in  the  collection  process  was  sought  in  order  to  attempt  to 
understand  the  limitations  involved  in  the  compensation  process  of  3D  motion  IMU  data.  This 
involved  sweeping  the  IMU  sensor  by  itself  from  side-to-side  and  placing  it  on  the  ground  at  the 
end  of  each  sweep  for  zero  velocity  updates.  To  implement  this,  a  series  of  fixed  locations  on  the 
test  platform  were  marked  by  right  angle  metal  braces  and  held  in  place  using  wooden  dowels 
(see  Figure  50).  The  locations  were  numbered  1  through  6  and  midpoint  to  midpoint  separations 
noted  so  that  the  relative  zero  velocity  update  positions  were  known.  The  IMU  was  then  simply 
lifted  and  hand-placed  at  a  series  of  locations,  with  the  data  being  collected  via  a  serial  cable 
onto  a  laptop  running  Gyro  View. 

The  IMU  was  always  oriented  the  same  way  on  the  left  side  (i.e.  locations  1,  3  &  5)  and  90 
degrees  from  this  on  the  right  side  (i.e.  locations  2,  4  &  6).  The  start  and  finish  IMU  positions 
were  always  at  location  1.  It  should  be  noted  that  the  IMU  placement  was  not  always  exact  - 
sometimes  it  landed  on  the  metal  brace  and  was  slightly  tilted,  and  sometimes  it  was  up  to  %  “ 
off  position.  To  average  these  out,  the  IMU  was  moved  back  and  forth  between  locations 
multiple  times,  and  the  measurements  repeated  several  times. 

6.3.  Processing  and  analysis 

The  zero  velocity  update  method  comprises  of  the  following  sequence  of  steps: 

•  Determination  of  the  angular  rate  bias  vector 

•  Determination  of  the  acceleration  bias  vector 

•  Determination  of  the  orientation  angles  of  the  IMU  body  relative  to  the  Earth’s  frame 

•  Integration  of  the  fully  compensated  acceleration  vector  in  the  Earth’s  frame 
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Figure  50:  Setup  for  3D  motion  of  just  the  IMU. 


The  qualifier  “fully”  implies  compensation  of  both  the  bias  vector  and  the  gravitational 
contribution,  over  all  time.  In  the  following  sections  we  will  trace  these  processing  steps  using 
the  3D  IMU  data  collected  as  described  in  the  discussion  of  Figure  50.  For  comparison,  we  will 
compare  the  results  with  performance  based  on  processing  static  IMU  data. 

6.3.1.  Determination  of  the  Angular  Rate  Bias  Vector 


The  determination  of  the  angular  rate  bias  vector  via  (20)  is  relatively  straightforward.  With 
sufficient  zero  velocity  updates  taken  periodically,  a  global  representation  of  b a(t)  can  be 
defined.  To  first  order,  this  is  well  represented  by  separate  linear  trends  for  each  of  the  three 
components,  i.e. 

MO 
MO 
MO 

The  Static  IMU 

Consider  a  50  minute  time  series  of  static  IMU  data.  Performing  a  linear  regression  on  each  of 
the  angular  rate  components  over  the  50  minute  segment,  yields 


=  m  t  +  c 

(Ox  CO 

-  m  t  +  cr 

a  ct 

=  m„  t  +  c„ 
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ma  =11.47 xl0“6  °/s2  and  ca  =  0.19880433  75 
mco  =0.60x10  6  °/s2  and  ^=0.11388729  75 
m  =  1.39x  10~6  °/ s2  and  c  =0.11041257  75 

C0Z  (Oz 

Over  the  entire  50  minutes  of  data,  the  compensated  angular  rate  components  give 

( ax(t)-bco  {t )}  =  -8.0xl0“14  7 5  and  a(no> )  =  0.300  °/s 
(cdy{t)-b0)  {t))  =  -5.1xWX5  °/s  and  a(na  )  =  0.332  75 
(®z(0-\(0)  =  -2.6xlO“16  75  a(na  )  =  0.327  °/ s 

where  (•}  denotes  averaging,  and  cr(»)  denotes  computing  the  standard  deviation. 

If  instead,  consecutive  1  second  windows  of  data  over  the  entire  50  minutes  of  data  are 
considered,  we  have 

(ax ( t )  -  b0>  (0)  =  5.85  x  10'7  ±  0.01 0465650  7 5  and  a(nco  )  =  0.301  +  0.020  7 5 
(cdy ( t )  - bm  (o)  =  -1 .89  X  10“7  ±  0.014537236  75  and  cr(nco  )  =  0.332 ±  0.023  7 5 
( S):{t)-ba  (t )}  =  -3.75 x  10'7  ± 0.014223592  75  rr(nfo )  =  0.327 ±  0.022  7  5 

It  should  be  noted  that  although  the  mean  of  the  average  values  still  approaches  zero,  local 
deviations  on  the  order  of  0.01  °/s  occur  as  a  result  of  the  linear  trend  assumptions  made. 

The  Unconstrained  Moving  IMU 

Consider  the  IMU  3D  motion  data  for  the  back  and  forth  sweeping  between  two  fixed  points 
defined  by  right  angle  metal  braces  secured  to  the  test  platform.  Zero  velocity  updates  for  a 
minimum  of  5  seconds  each  were  conducted  at  the  fixed  points.  The  time  in  motion  between  the 
zero  velocity  updates  typically  lasted  under  3  seconds. 

Performing  a  linear  regression  on  each  of  the  angular  rate  components  (using  the  data  over  all 
zero  velocity  update  segments),  yields 

ma  =  -5.15  x  10‘6  7 52  and  c ^  =  0.16389530  7 5 
m  =  -123.13 x  10~6  7 s1  and  c  =0.24644612  75 

COy  COy 

ma  =  67.94  x  10“6  7 52  and  ca  =  0.08594225  7 5 

Comparing  these  regression  coefficients  to  the  ones  obtained  in  the  previous  example,  it  is  clear 
that  both  the  slope  and  initial  bias  coefficient  values  can  vary  significantly.  This  reinforces  the 
importance  of  taking  zero  velocity  updates,  if  for  no  other  reason  than  to  detennine  the 
appropriate  b a(t)  for  the  data  in  question. 
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Figure  5 1  below  shows  the  mean  (diamond  symbols)  and  standard  deviation  (error  bars)  of  the 
angular  rate  components  at  each  of  the  zero  velocity  updates,  with  the  derived  linear  regression 
curves  over-plotted  in  red.  The  mean  values  have  been  placed  at  the  temporal  midpoint  of  each 
zero  velocity  update  duration. 
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Figure  51:  The  angular  rate  bias  vector  components  plotted  (in  red]  with  the  mean  and  standard  deviation 
for  each  of  the  zero  velocity  updates. 

Subtracting  each  linear  regression  curve  from  the  respective  angular  rate  component,  gives,  over 
all  zero  velocity  update  segments 


(S  (0  -  ba  (0)  =  -2. 1 1 X 1 0“7  +  0.004490364  7  5  and  a(noj  )  =  0.305  ±  0.01 3  °/s 
(coy(t)  -  ba  (0)  =  -1 .22103  x  10“4  ±  0.005412638  7 s  and  a(na> )  =  0.326  ±  0.01 1  °/s 

(S  (t)  -  bo>  (0)  =  -1 .252090  x  1 0“3  ±  0.0047 1 41 05  7  5  and  cr(nco  )  =  0.322  ±  0.01 1  7  j 

which  is,  more  or  less,  in  line  with  what  was  obtained  for  the  compensated  angular  rate 
components  in  the  previous  example  of  the  static  IMU. 


6.3.2.  Determination  of  the  Acceleration  Bias  Vector 
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The  determination  of  the  acceleration  bias  vector  via  (21)  is  more  ambiguous.  During  any 
particular  zero  velocity  update,  the  term  Re^(t)g  will  contribute  an  orientation-dependent 
constant  bias  vector  to  the  measured  acceleration  vector.  Without  an  independent  measure  of  this 
contribution,  the  best  we  can  do  to  determine  the  acceleration  bias  vector  is  to  use  only  zero 
velocity  update  data  where  the  IMU  orientation  was  known  to  be  identical.  In  so  doing,  (16)  can 
be  recast  as 

(27)  a(0  =  bl(0  +  nfl(0 

where  the  global  representation  of  b^t)  now  has  the  orientation-dependent  constant  vector 
absorbed  into  it.  As  with  the  angular  rate  bias  vector,  to  first  order,  b',(0  is  well  represented  by 
separate  linear  trends  for  each  of  the  three  components,  i.e. 

b\_,  (t)  =  >n„  t  +  ca 

b\y  (0  =  ma,  t  + 

b' ai  if)  =  mj  +  c3i 

where  now  cfl  =  R,,  >A(0)g  +  b„(0) .  Again,  since  we  do  not  have  independent  confirmation  of 
R,,_>A(0)g  ,  it  is  currently  assumed  that  b  ,(0)  =  0  . 

The  Static  IMU 

Consider  the  same  50  minute  time  series  of  static  IMU  data  as  before.  Noting  that  the  orientation 
of  the  IMU  body  remains  unchanged  throughout,  a  linear  regression  on  each  of  the  acceleration 
components  over  the  entire  50  minutes  can  be  performed,  yielding 

ma  =9.9x10  ^g/s  and  ca  =  0.000171506g 
ma  =  -60.2  x  10  s  gl s  and  ca  =  0.020438798g 
ma  =  0.5  x  10“8g/j  and  ca  =  0.998410810g 

Over  the  entire  50  minutes  of  data,  the  compensated  acceleration  components  give 

(ax(t)  -  b'a  (0)  =  3.8  x  10~16  g  and  cr(na  )  =  0.001 159g 
(ay(t)-b'a  (t)^  =  8.6x1 0~15g  and  cr(»a  )  =0.001390g 
(a_{t)  -  b'a  (tfj  =  1.6  x  10-13  g  and  <j(na  )  =  0.001 174g 

If  instead,  consecutive  1  second  windows  of  data  over  the  entire  50  minutes  of  data  are 
considered,  we  have 
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(ax  (t)  -  b'a  (tj)  =  3 .0  x  1 0“9  ±  0.000 1 29 1 83  5  g  and  a(na  )  =  (1 1 . 52  ±  0.95)  x  1  O'4  g 
(dy(t)-b'a  (O)  =  1.5x1 0“9  +0.0001 280888 g  and  cj(na  )  =  (13.85  +  1.00)  xl O'4  g 
(a,(t)  -  b[  (0)  =  -3.0  x  10‘9  ±  0.0001 145845g  and  o{na  )  =  (11.69  +  0.94)  x  10“4  g 

Again,  although  the  mean  of  the  average  values  still  approaches  zero,  local  deviations  on  the 
order  of  0.1  mg  occur  as  a  result  of  the  linear  trend  assumptions  made. 

The  Unconstrained  Moving  IMU 

Consider  the  same  IMU  3D  motion  data  as  before.  Performing  a  linear  regression  on  each  of  the 
acceleration  components,  using  only  the  initial  and  final  zero  velocity  update  data  (the  IMU 
orientation  was  known  to  be  identical  during  these  updates),  yields 

ma  =  -13.75  x  10~7  g/s  and  ca  =  0.000434544 g 
ma  =  3.82  x  10~7  g/s  and  ca  =  0.001 184215g 
mu  =-3.23x10 -1  g/s  and  ca  =  0.998566701g 

Subtracting  each  linear  regression  curve  from  the  respective  acceleration  component,  gives,  over 
the  initial  and  final  zero  velocity  update  segments 

(ax(t)  -  b[  (; t )}  =  (4.88  ±  1 1.82)  x  10“7  g  and  a(na  )  =  (9.14  ±  0.64)  x  10~4  g 
(ay(t)  -  b'a  (0)  =  (36.87  ±  89.28)  x  10‘7  g  and  cr(na  )  =  (14.68  ±  0.52)  x  10“4  g 
(az(t)  -  b'a  (; t ))  =  (1 1.65  ±  28.23)  x  10  1  g  and  o(na  )  =  (11.35  +  2.03)  x  10“4  g 

The  fact  that  the  deviations  about  the  mean  are  much  smaller  than  in  the  previous  example  is  a 
direct  consequence  of  having  an  inadequate  statistical  sampling.  To  be  sure,  if  all  the  zero 
velocity  update  data  where  the  IMU  orientation  was  known  to  be  identical  (not  just  the  initial 
and  final  segment  data)  is  instead  used,  we  have  for  the  new  linear  regression  coefficients 

ma  =-16.04x10 ~7  g/s  and  ca  =  0.0004638 14 g 
ma  =  7.16  x  10  7  g/s  and  ca  =  0.001 1 19055 g 
ma  =  -3.56  x  10~7  g/s  and  ca  =  0.998572324 g 

Subtracting  each  linear  regression  curve  from  the  respective  acceleration  component,  gives,  over 
the  same  zero  velocity  update  data  used  to  perfonn  the  regressions 
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(ax{t)  -  b'a  (0)  =  (0.86  ±  7.86)  x  10“5  g  and  a(na  )  =  (8.95  ±  0.42)  x  10 ^ 

(dy(t)-b'a  (I)}  =  (-5.23  ±10.46)  x  10'5  g  and  a(na  )  =  (15.1 1  ±  0.54)  x  10“4  g 
(a,(t)  -  b[  (*))  =  (0.55  ±  2.25)  x  10‘5  g  and  a(na  )  =  (10.78  ±  1 .14)  x  10‘4  g 

which  is  again,  more  or  less,  in  line  with  what  was  obtained  for  the  compensated  acceleration 
components  in  the  previous  example  of  the  static  IMU. 

6.3.3.  Determination  of  the  Orientation  Angles 

As  noted  earlier  (19),  the  attitude  of  the  IMU  body  frame  relative  to  the  Earth  frame  is 
determined  from  the  compensated  angular  rates  co,  (t) ,  cov{t)  and  coAt) ,  via  the  following 
coupled  differential  equations 

<f>(t)  1  sin  (j){t)  tan  0{t)  cos  0(f)  tan  6*(t)  cox(t) 

—  0(t)  =  0  cos  0(f)  -sin0(t)  cov(t) 

dt\y(t)\  o  Sin</>{t)  C0S</>{t)  [^(0 

cos6*(t)  cos6*(t) 

where  y/(t),  0(t)  and  (/At)  are  the  yaw,  pitch  and  roll  angles,  respectively.  During  zero  velocity 
updates,  an  additional  independent  measure  of  both  the  pitch  and  roll  angles  is  possible  by 
making  use  of  the  de-biased  acceleration  data.  Specifically,  averaging  (21)  over  the  extent  of  the 
7-th  zero  velocity  update  that  starts  at  time  r, ,  gives 

(28)  (a(t)-bfl(0)  =  R^(fi)g 

where,  explicitly  in  terms  of  the  Euler  angles, 

-  sin  6t 

(29)  Re->A(E)g=  cos 0z  sin 0.  g 

cos  0T  COS  0r 

and  the  constant  angles  0T  and  0r  specify,  respectively,  the  pitch  and  roll  of  the  IMU  body 
frame  relative  to  the  Earth  frame  during  the  7-th  zero  velocity  update. 

From  (28)  and  (29), 

(ax(t)-ba  (0)  I"  -sin^r 

(30)  (fi  (0-^,  (0)  =  cos^  sin0r  g 
(dz(t)-ba(t))  [cos^(cos0ri_ 

allowing  one  to  determine  the  angles  via 
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(31) 


(flx(0-\(0)  1 


0 


=  sin 


V 


g 


and 

(32) 


(j)T  =  sin- 


f  (ay(t)-ba(tf^ 


g  cos  6t 


=  cos 


£(0~M0) 

g  cos  0r 


It  is  implied  in  (31)  and  (32)  that  the  averaging  takes  place  over  the  duration  of  the  z'-th  zero 
velocity  update. 

The  Static  IMU 

Again,  consider  the  same  50  minute  time  series  of  static  IMU  data  as  before.  Using  the 
compensated  angular  rates,  and  assuming  zero  initial  angles,  the  pitch  and  roll  angles  computed 
via  (19)  are  shown  in  Figure  52  below. 


Figure  52:  The  pitch  and  roll  angles  determined  via  (19),  for  the  case  of  the  IMU  sitting  still  on  a  table  top 
for  50  minutes.  Compensated  angular  rates  have  been  used  and  initial  angles  of  zero  assumed. 


Now  assume  that  we  have  zero  velocity  updates  of  5s  duration,  followed  by  3s  of  unknown 
motion,  for  the  entire  50  minutes.  The  zero  velocity  updates  provide  an  additional  independent 
measure  of  both  pitch  and  roll  angles,  via  (31)  and  321).  Since  the  computed  angles  over  all  the 
zero  velocity  updates  vary  by  less  than  0.01°,  these  angles  serve  well  to  constrain  inevitable 
angular  errors  that  evolve  via  the  usage  of  (19). 

Figure  53  shows  the  result  of  forcing  the  pitch  and  roll  angles  over  the  zero  velocity  updates 
(shown  in  red)  to  be  those  computed  via  (31)  and  (32),  while  using  (19)  to  determine  the  angles 
between  the  zero  velocity  updates  (shown  in  black).  The  additional  effect  of  removing  linear 
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trends  between  the  zero  velocity  updates,  so  that  angles  at  the  end  of  these  segments  agree  with 
angles  at  the  subsequent  zero  velocity  update,  is  shown  in  Figure  54.  The  bottom  two  panels  in 
both  Figures  53  and  54  show  close  ups  of  the  first  200  seconds  of  the  data  in  the  top  two  panels. 

The  Unconstrained  Moving  IMU 

Again,  consider  the  same  IMU  3D  motion  data  as  before.  Using  the  compensated  angular  rates, 
and  assuming  zero  initial  angles,  the  pitch  and  roll  angles  computed  via  (19)  are  shown  in  Figure 
55  below.  The  data  at  the  two  fixed  points  (about  which  the  sweeping  takes  place)  are  explicitly 
represented  in  red,  for  one  point,  and  green,  for  the  other.  Apart  from  an  obvious  anomalous 
orientation  on  the  second  return  to  the  first  point,  gradual  temporal  trends  are  observed  for  the 
angles  at  the  fixed  points  -  most  notably  for  the  pitch  angle  at  the  first  fixed  point.  The  trends, 
however,  are  not  well  represented  by  simple  low  order  polynomials. 

Figure  56  shows  the  result  of  forcing  the  pitch  and  roll  angles  over  the  zero  velocity  updates 
(taken  at  the  fixed  points  mentioned  above)  to  be  those  computed  via  (31)  and  (32),  while  using 
(19)  to  determine  the  angles  between  the  zero  velocity  updates.  The  additional  effect  of 
removing  linear  trends  between  the  zero  velocity  updates,  so  that  angles  at  the  end  of  these 
segments  agree  with  angles  at  the  subsequent  zero  velocity  update,  is  shown  in  Figure  57. 

By  using  compensated  pitch  and  roll  angles,  along  with  compensated  angular  rates,  a  favorable 
estimate  of  the  yaw  angle  via  (19)  is  now  possible.  The  top  panel  of  Figure  58  shows  the  yaw 
angle  obtained  by  simply  assuming  an  initial  angle  of  -45°  in  (6).  Since  the  IMU  body  is  known 
to  eventually  return  to  its  initial  position  and  orientation,  a  final  yaw  angle  of  -45°  is  also 
assumed.  The  bottom  panel  shows  the  result  of  subtracting  a  linear  trend  from  the  yaw  angle  in 
the  top  panel,  so  that  the  final  computed  angle  agrees  with  the  assumed  final  angle  of  -45°. 
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Figure  53:  The  pitch  and  roll  angles  determined  via  (31)  and  (32),  respectively,  using  the  de-biased 
accelerations  over  assumed  5  second  zero  velocity  updates  |in  red],  and  via  (19)  using  the  compensated 
angular  rates  during  assumed  3  second  segments  of  unknown  motion  between  updates  [in  black],  for  the  case 
of  the  IMU  sitting  still  on  a  table  top  for  50  minutes.  The  bottom  two  panels  are  close  ups  of  the  first  200 
seconds. 
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Figure  54:  The  pitch  and  roll  angles  of  Figure  53,  but  with  linear  trends  subtracted  from  the  assumed 
segments  of  unknown  motion  between  updates.  Any  discontinuities  that  existed  between  the  angles  at  the  end 
of  these  segments  and  those  at  the  subsequent  zero  velocity  updates  have  been  removed.  As  before,  the 
bottom  two  panels  are  close  ups  of  the  first  200  seconds. 
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Figure  55:  The  pitch  and  roll  angles  determined  via  (19),  for  the  case  of  back  and  forth  3D  sweeping 
motion  of  the  IMU  between  two  well-defined  fixed  points  (represented  by  the  red  and  green  data]. 
Compensated  angular  rates  have  been  used  and  initial  angles  of  zero  assumed. 


Figure  56:  The  pitch  and  roll  angles  determined  via  (31)  and  (32),  respectively,  using  the  de-biased 
accelerations  during  the  zero  velocity  updates  [in  red  and  green],  and  via  (19)  using  the  compensated  angular 
rates  between  updates  (in  black],  for  the  case  of  back  and  forth  3D  sweeping  motion  of  the  IMU  between  two 
well-defined  fixed  zero  velocity  update  points. 
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Figure  57:  The  pitch  and  roll  angles  of  Figure  56,  but  with  linear  trends  subtracted  between  zero  velocity 
updates  so  as  to  remove  any  existing  discontinuities. 
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Figure  58:  The  yaw  angle  determined  via  (19)  using  the  compensated  angular  rates  and  the  compensated 
pitch  and  roll  angles  of  Figure  57.  The  yaw  angles  during  the  initial  and  final  zero  velocity  updates  were 
assumed  to  be  -45°.  The  bottom  panel  shows  the  compensated  yaw  angle,  where  a  linear  trend  has  been 
subtracted  between  the  end  of  the  initial,  and  the  beginning  of  the  final,  zero  velocity  updates. 


6.3.4.  Integration  of  the  Fully  Compensated  Acceleration  Vector 
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Having  determined  the  compensated  pitch  and  roll  angles,  the  gravitational  contribution  to  the 
acceleration  vector,  RMj(f)g ,  is  now  known  over  all  t.  With  this,  an  estimate  of  the  true  (or 
fully  compensated)  acceleration  vector  a (t)  is  achievable  via  (16). 

Recall  that  the  vector  a (t)  of  (16)  is  defined  in  the  IMU  body  frame.  Since  the  goal  is  to 
determine  the  position  of  the  IMU  in  the  Earth  frame,  it  is  necessary  to  transform  the  vector  a (t) 
from  the  IMU  body  frame  to  the  Earth  frame  before  carrying  out  the  successive  integration 
operations  to  obtain  the  velocity,  followed  by  the  position,  vectors. 

The  Static  IMU 

Consider  the  same  static  IMU  data  as  before  and,  again,  assume  that  we  have  zero  velocity 
updates  of  5s  duration,  followed  by  3s  of  unknown  motion,  for  the  entire  50  minutes.  The 
compensated  pitch  and  roll  angles  have  already  been  determined  for  this  case  (refer  to  Figure 
54).  The  gravitational  contribution  to  the  acceleration  vector,  for  all  t,  is  therefore  known  and  the 
full  compensation  of  the  acceleration  vector  is  now  possible  via  (16).  Figure  59  below  shows 
that  when  this  is  evaluated,  however,  a  residual  bias  in  the  z-component  still  persists,  requiring  a 
revision  in  our  previous  assumption  of  ba(0)  =  0.  We  shall  see  later  that  a  more  thorough 
revision  of  this  assumption  is  needed  in  order  to  attain  the  cm  level  accuracy  in  positioning  that 
this  project  ultimately  seeks. 

Subtracting  the  residual  bias  from  the  z-component  yields  the  fully  compensated  acceleration 
vector.  Transforming  this  vector  from  the  IMU  body  frame  to  the  Earth  frame  and  integrating, 
results  in  the  velocity  vector  shown  in  Figure  60  below.  Since,  by  definition,  the  velocity  is  zero 
during  zero  velocity  updates,  the  integration  process  is  only  carried  out  between  updates,  during 
segments  of  assumed  unknown  motion.  The  additional  effect  of  subtracting  linear  trends  between 
the  zero  velocity  updates  is  shown  in  Figure  61  below.  This  is  done  so  as  to  force  agreement 
between  what  the  integration  tells  us  the  velocity  ought  to  be,  and  what  we  know  the  velocity  to 
be,  at  the  subsequent  updates. 

Finally,  integrating  the  compensated  velocity  vector  of  Figure  61  provides  an  estimate  of  the 
position  vector.  Figures  62  and  63  show  the  outcome  of  this  integration  process  with  three 
separate  assumptions  used.  If  no  assumption  other  than  an  initial  position  of  zero  is  made,  the 
results  are  the  black  curves  of  Figures  62  and  63.  In  order  to  reduce  the  positional  errors  that 
inevitably  accumulate  over  time,  assumptions  based  on  revisiting  known  positions  -  i.e.  zero 
position  updates,  in  this  case  -  are  examined.  The  green  curves  represent  the  result  of  assuming  a 
final  zero  position  update,  while  the  red  curves  represent  the  result  of  assuming  zero  position 
updates  at  every  25th  zero  velocity  update.  It  is  clear  from  Figure  63  that  with  sufficient  zero 
position  updates,  the  errors  can  be  reduced  to  the  desired  cm  level  for  the  entire  50  minute 
duration. 
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Figure  59:  The  mean  and  standard  deviation  of  the  compensated  acceleration  components  over  each  of  the 
assumed  zero  velocity  updates,  for  the  case  of  the  IMU  sitting  still  on  a  table  top  for  50  minutes.  Due  to  the 
proximity  and  shear  quantity  of  updates,  the  diamonds  representing  the  means  appear  to  merge  into 
continuous  thick  dark  curves,  while  the  error  bars  produce  the  lighter  fabric  that  envelope  the  mean  curves. 
The  red  curves  represent  residual  linear  trends. 
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Figure  60:  The  velocity  components  obtained  by  integrating  the  corresponding  fully  compensated 
acceleration  components  between  the  assumed  zero  velocity  updates,  for  the  case  of  the  IMU  sitting  still  on  a 
table  top  for  50  minutes.  By  definition,  the  velocity  is  zero  during  the  zero  velocity  updates.  The  fully 
compensated  acceleration  components  used  here  have  had  the  residual  trends,  shown  in  Figure  59, 
subtracted  from  each  component,  and  have  also  been  transformed  from  the  IMU  body  frame  to  the  Earth 
frame. 
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Figure  62:  The  estimated  position  components  in  the  Earth  frame  obtained  by  integrating  the 
corresponding  compensated  velocity  components  of  Figure  61.  The  black  curves  represent  direct 
integrations,  with  initial  positions  of  zero  being  the  only  assumptions  made;  the  green  curves  are  the  result  of 
forcing  the  final  positions  to  be  zero  as  well,  with  linear  corrections  applied  between  both  ends;  and  the  red 
curves  are  the  result  of  forcing  every  25th  update  (200  seconds)  to  be  zero,  with  linear  corrections  applied 
between  each  of  these. 
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Figure  63:  The  estimated  X  and  Y  positions  of  Figure  62  plotted  against  one  another,  showing  the 
deviations  from  the  true  (0,0)  position.  The  blue  circle  shown  has  a  radius  of  2  cm  about  the  true  position. 
The  same  color  coding  used  in  Figure  62  is  retained  here. 


The  Unconstrained  Moving  IMU 

Again,  consider  the  same  IMU  3D  motion  data  as  before.  Analogous  to  the  static  IMU  case,  the 
evaluation  of  the  full  compensation  of  the  acceleration  vector  via  (16)  reveals  the  presence  of  a 
residual  bias  in  the  z-component.  As  shown  in  Figure  64,  the  magnitude  of  this  bias  is  similar  to 
that  observed  in  Figure  59,  but  not  identical. 

The  source  of  the  residual  bias  is  evident.  Recall  the  assumption  made  (i.e.  ba(0)  =  0  )  following 
the  introduction  of  (27)  so  that  the  linear  regressions  on  zero  velocity  update  acceleration  data 
with  identical  IMU  orientation  resulted  in  ca  =  Re  >6(())g .  Clearly  from  the  examples  that 

immediately  follow  this,  |ca|  ^  g ,  even  though  the  physics  require  it.  A  simple  way  to  resolve 
this  inconsistency  is  by  requiring  that  ba  (0)  ^  0  such  that  |ca  -  ba  (0)|  =  g  .  However,  while  this 
simple  method  should  directly  account  for  the  residual  bias  discussed  above,  it  does  not. 
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Figure  64:  The  mean  and  standard  deviation  of  the  compensated  acceleration  components  over  the  zero 
velocity  updates,  for  the  case  of  back  and  forth  3D  sweeping  motion  of  the  IMU  between  two  well-defined 
fixed  points.  The  zero  velocity  updates  take  place  at  the  two  fixed  points 


The  velocity  vector  obtained  by  integrating  the  fully  compensated  acceleration  vector  (i.e.  the 
compensated  acceleration  vector  obtained  via  (16)  minus  any  residual  biases)  is  shown  in  Figure 
65.  Again,  since  the  velocity  is,  by  definition,  zero  during  zero  velocity  updates,  the  integration 
process  is  only  carried  out  between  updates.  As  before,  the  results  of  the  integration  process 
between  updates  will  diverge  from  what  they  should  be  by  the  time  the  subsequent  updates  are 
reached.  The  magnitudes  of  these  divergences  differ  from  the  static  IMU  case  by  about  a  factor 
of  ten  (c.f.  Figures  60  and  65). 
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Figure  65:  The  velocity  components  obtained  by  integrating  the  corresponding  fully  compensated 
acceleration  components  between  the  zero  velocity  updates,  for  the  case  of  back  and  forth  3D  sweeping 
motion  of  the  IMU  between  two  well-defined  fixed  points.  By  definition,  the  velocity  is  zero  during  the  zero 
velocity  updates.  The  fully  compensated  acceleration  components  used  here  have  had  the  residual  trends, 
shown  in  Figure  64,  subtracted  from  each  component,  and  have  also  been  transformed  from  the  IMU  body 
frame  to  the  Earth  frame. 

The  compensated  velocity  vector  components  are  shown  in  Figure  66.  As  before,  the 
compensation  is  obtained  by  removing  linear  trends  between  the  zero  velocity  updates.  This  step 
ensures  that  the  integrated  final  velocities  between  updates  agree  with  the  subsequent  update 
velocities.  Integration  of  these  compensated  velocities  result  in  the  estimated  positions 
represented  in  Figures  67  and  68.  Similar  to  the  static  IMU  case,  three  different  outcomes  of  the 
estimated  positions  are  considered:  one  where  no  assumptions  other  than  an  initial  position  of 
zero  is  made  (black  curves);  one  where  the  additional  assumption  that  we  return  to  the  initial 
position  at  the  end  is  made  (green  curves);  and  one  where  the  assumption  is  made  that  we  return 
to  the  initial  position  at  every  other  zero  velocity  update  (red  curves). 
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Figure  66:  The  velocity  components  of  Figure  65,  but  with  linear  trends  subtracted  between  zero  velocity 
updates  so  as  to  remove  any  existing  discontinuities. 

The  fixed  end  points  of  the  3D  sweeping  motion  are  represented  in  Figure  68  by  blue  circles  of  2 
cm  radius.  It  is  clear  from  the  top  panel  of  Figure  68  that  systematic  trends  in  both  the  x  and  y 
positions  exist.  However,  even  after  applying  the  evermore  restrictive  zero  position  update 
assumptions  (i.e.  lower  two  panels),  the  project  goal  of  cm  level  positional  accuracy  is  not  being 
met.  Even  with  position  updates  between  each  back  and  forth  sweep  (red  curves)  the  location 
accuracy  drifts  off  by  several  centimeters  over  the  course  of  the  sweeping  motion. 
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Figure  67:  The  estimated  position  components  in  the  Earth  frame  obtained  by  integrating  the 
corresponding  compensated  velocity  components  of  Figure  61.  The  black  curves  represent  direct 
integrations,  with  initial  positions  of  zero  being  the  only  assumptions  made;  the  green  curves  are  the  result  of 
forcing  the  final  positions  to  be  zero  as  well,  with  linear  corrections  applied  between  both  ends;  and  the  red 
curves  are  the  result  of  forcing  every  other  update  position  to  be  zero,  with  linear  corrections  applied 
between  each  of  these. 
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Figure  68:  The  estimated  X  and  Y  positions  of  Figure  67  plotted  against  one  another,  showing  the 
deviations  from  the  two  fixed  points  (blue  circles)  about  which  the  sweeping  takes  place.  The  blue  circles 
shown  have  radii  of  2  cm.  The  same  color  coding  used  in  Figure  67  is  retained  here. 
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7.  SAINT  Tests 


Although  we  had  good  results  in  determining  sensor  location  from  the  Crossbow  IMU  data  for 
2D  motion,  the  results  for  3D  motion  were  discouraging.  We  believe  that  this  is  because  of 
unacceptable  levels  of  noise  and  drift  in  the  Crossbow  angle  rate  measurements.  Accurate 
determination  of  the  IMU's  orientation  is  vital  to  getting  accurate  position  information  because 
even  a  very  small  error  in  angle  can  introduce  erroneous  gravity  contributions  to  the  X  and  Y 
accelerations,  which  in  turn  integrate  up  to  large  position  errors.  A  higher  quality  Honeywell 
HG1900  IMU  was  integrated  with  the  EM61-HH  by  ENSCO  under  ESTCP  project  MM-0604.  It 
is  referred  to  as  the  Small-Area  Inertial  Navigation  Tracking  or  SAINT  system.  We  collaborated 
with  MM-0604  to  evaluate  whether  or  not  the  SAINT  could  provide  sufficiently  accurate  sensor 
position  information  to  support  reliable  UXO/clutter  discrimination.  Results  are  reported  in  [9]. 

7.1.  The  SAINT  system 

The  SAINT  system  is  shown  in  Figure  69.  It  comprises  a  Honeywell  HG1900  IMU,  Leica  digital 
magnetic  compass  and  embedded  CPU  to  record  the  data  to  compact  flash  card.  The  EM61-HH 
data  is  collected  by  the  standard  Geonics  Allegro  logger.  Calculation  of  sensor  location  from  the 
IMU  is  done  in  post-processing. 


Allegro  Bracket 

SAINT  Enclosure 
with  IMU 


Tripod 
EM61-HH  Coil 


Figure  69:  SAINT  System  combining  EM61-HH  and  Honeywell  HG1900  IMU. 

The  HG1900  is  a  more  expensive  unit  (~$  1  OK  vs.  ~$4K),  but  its  angular  rate  random  walk  spec 
(0.1  deg/hr12)  is  substantially  better  than  that  of  the  Crossbow  IMU400CC-100  (2.25  deg/hr12). 
It  also  has  a  significantly  higher  bandwidth  than  the  Crossbow  unit  (600  Hz  vs.  135  Hz).  Figure 
70  shows  simultaneous  angle  rate  data  from  the  Crossbow  unit  and  the  Honeywell  unit.  They 
were  mounted  side-by-side  on  the  shaft  of  the  EM61-HH.  The  data  in  Figure  70  were  collected 
as  the  sensor  was  being  swept  back  and  forth.  Several  things  are  apparent  in  these  data:  the 
Honeywell  data  appears  smoother  (higher  bandwidth)  and  less  noisey  than  the  Crossbow  data, 
and  there  appears  to  be  a  bit  of  drift  with  the  Crossbow  unit,  especially  in  the  Y-axis  angle  rate. 
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Figure  70:  Simultaneous  angle  rate  data  from  Crossbow  and  Honeywell  IMUs. 


7.2.  Test  Setup 

As  part  of  ESTCP  project  MM-0604  the  SAINT  IMU/EM61-HH  was  used  to  interrogate  a 
selection  of  the  Calibration  Grid  buried  inert  ordnance  at  the  Aberdeen  Proving  Ground 
Standardized  Test  Site.  The  EM61-HH  sensor  was  swept  back  and  forth  over  known  object 
locations  and  positioning  for  this  3-D  sensor  trajectory  was  measured  using  the  Honeywell  unit. 
We  set  up  the  Arcsecond  system  to  provide  ground  truth  on  a  subset  of  this  data.  For  the  data 
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where  the  EM61-HH  sensor  was  simultaneously  tracked  with  the  Arcsecond  system  and  the 
IMU,  we  found  that  typically  the  agreement  between  the  IMU-derived  sensor  positions  and  the 
Arcsecond  ground  truth  was  good  to  within  1-3  cm. 

The  Aberdeen  Calibration  Grid  is  set  up  into  rows  A  through  K  and  columns  1  through  15 
increasing  from  the  NE  corner  of  the  grid.  The  smaller  buried  items  are  separated  by  2  meters. 
For  this  test,  a  selection  of  items  at  shallower  depths  was  interrogated,  because  of  the  limitations 
on  the  EM61-HH’s  detection  of  deep  objects  and  the  overly  challenging  depths  used  on  the 
Calibration  Grid.  The  items  covered  were  20mm,  40mm  grenade,  57mm,  60mm  mortar,  2.75 
inch  rocket  warhead,  81mm  mortar,  105mm,  155mm,  and  a  shot-put.  For  the  most  part,  three 
examples  of  each  item  were  interrogated  using  only  the  ENSCO  Honeywell  system.  Each 
interrogation  over  a  given  object  was  repeated  at  least  three  times,  resulting  in  three  independent 
data  sets  over  each  object.  A  20mm  at  the  grid  location  A3  was  measured  repeatedly  with  both 
the  IMU  and  the  Arcsecond  system.  In  this  case,  five  measurements  were  taken  in  the  morning 
and  then  five  more  in  the  afternoon. 

The  Arcsecond  system  that  we  used  included  a  set  of  four  fixed  transmitter  stations  and  the  four 
receivers  attached  to  the  EM61-HH.  After  a  calibration  set  up,  the  system  collects  data  to  track 
the  EM61-HH  coil  head  in  3-D  along  with  orientation  angles  of  yaw,  pitch,  and  roll.  The  actual 
full  trajectory  at  a  20  Hz  data  rate  is  only  available  in  post-processing  of  the  receiver 
measurements.  The  current  post-processing  only  provides  the  trajectory  in  a  local  Arcsecond 
frame  of  reference  based  on  the  location  of  the  fixed  transmitters.  In  order  to  map  the  Arcsecond 
data  to  either  the  local  APG  grid  or  to  UTM  coordinates  a  surveyed  point  in  the  northeast  corner 
of  the  site  and  a  set  of  known  grid  points  were  measured  with  the  Arcsecond  system. 

The  ENSCO/Honeywell  inertial  system  integrates  measurements  from  the  HG1900  IMU.  The 
integration  starts  and  stops  from  a  given  “zero  point”.  All  data  is  relative  to  this  zero  point.  The 
direction  of  the  ENSCO  local  coordinates  is  determined  by  the  Leica  electronic  compass 
included  in  the  system.  In  order  to  interrogate  all  of  the  cells  across  the  site,  four  different  zero 
points  were  used.  To  translate  from  ENSCO’s  local  coordinates  to  Arcsecond’s,  these  zero 
points  were  measured  with  the  Arcsecond.  The  IMU  raw  data  is  collected  at  600  Hz.  After  post¬ 
processing  this  data  to  an  actual  trajectory  (x,  y,  z,  yaw,  pitch,  and  roll),  ENSCO  down  sampled 
the  position  results  to  20  Hz.  The  ENSCO  data  was  provided  in  a  positive  z  pointing  down 
coordinate  system  and  had  to  be  converted  to  match  the  upward  z  pointing  system  used  by  the 
Arcsecond  and  our  EM  fitting  algorithms. 
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Figure  71:  EM61-HH  coil  head  trajectory  over  A3  cell  in  ArcSecond  coordinate  system.  IMU-derived 
trajectory  is  plotted  in  black  and  ArcSecond  trajectory  in  red. 

7.3.  Trajectory  Comparisons 

To  compare  the  two  measured  trajectories,  the  ENSCO  data  was  shifted  in  time,  shifted  in  space, 
and  rotated  to  best  match  the  Arcsecond  trajectory.  The  time  shift  is  due  to  the  two  different  data 
acquisition  clocks  used  by  the  systems.  The  spatial  shift  was  based  on  the  ENSCO  zero  point 
measured  by  the  Arcsecond.  For  the  ten  common  data  sets,  roughly  the  same  rotation  angle 
provided  a  good  match  between  the  two  sets  of  data.  The  average  angle  was  41.2  degrees 
counter-clockwise  to  match  the  ENSCO  data  to  the  Arcsecond.  The  angles  ranged  from  40.1  to 
42.05  with  a  standard  deviation  of  0.6.  Figure  71  plots  a  comparison  of  the  two  trajectories  in  the 
X,  Y  plane  for  one  of  the  joint  data  collections  over  the  A3  cell  location.  The  black  curve  is  from 
the  ENSCO  system  and  the  red  is  from  the  Arcsecond.  The  red  plus  signs  are  measured  grid  and 
zero  point  locations. 

To  match  the  UTM  coordinate  system,  a  variety  of  local  grid  points  was  measured  by  the 
Arcsecond.  Based  on  these,  the  Arcsecond  coordinates  must  be  rotated  35.9  degrees  clockwise  to 
align  with  true  north.  From  this,  the  ENSCO  coordinates  needed  to  be  rotated  5.3  degrees 
counter-clockwise  to  match  UTM. 

Using  the  average  rotation  angle  and  zero  point,  Figure  72  plots  the  ENSCO-only  trajectories 
from  the  same  zero  point.  This  time  the  data  is  in  the  local  grid  coordinates  relative  to  the 
northeast  corner  of  the  A1  cell.  The  red  curve  is  from  a  common  Arcsecond/ENSCO  trajectory. 
The  black  curves  are  all  from  ENSCO  only  data  collections.  The  red  plus  signs  are  fixed  grid 
points  measured  by  the  Arcsecond  system.  For  these  ENSCO  only  trajectories,  there  appears  to 
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be  some  further  rotation  and  possibly  shift  required  to  get  them  to  match  the  actual  measured  cell 
location.  Figure  73  plots  similar  results  from  two  other  zero  points. 

Figures  74  and  75  graph  comparisons  of  joint  ArcSecond/ENSCO  trajectories  in  more  detail.  In 
74(a)  and  75(a),  the  ArcSecond  trajectory  is  plotted  in  red  and  ENSCO  is  in  black.  The  other 
graphs  plot  (b)  the  horizontal  and  (c)  vertical  distances  between  the  two  trajectories  as  a  function 
of  time.  In  Figure  74  there  is  very  close  agreement  and  in  Figure  75  there  is  not.  The  Arcsecond 
system  typically  jitters  about  the  true  trajectory,  and  occasionally  spikes  way  off  when  one  or 
more  receivers  are  blocked.  However,  its  absolute  positioning  over  time  tends  to  be  very  good. 
The  ENSCO  data  plots  a  very  smooth  trajectory,  but  appears  to  drift  off  from  the  Arcsecond 
trajectory.  Out  of  the  ten  joint  trials,  eight  of  them  drifted  off  on  the  order  of  one  to  three 
centimeters  in  the  horizontal  similar  to  the  example  in  Figure  74.  For  two  of  the  trials,  the 
ENSCO  trajectory  drifted  way  off,  and  the  extent  of  the  drift  varied  over  time  as  in  Figure  75. 
Comparing  the  orientation  angles  between  the  two  showed  similar  results  with  typical 
differences  in  orientation  being  on  the  order  of  only  several  degrees  or  less.  Statistics  for  the 
differences  in  position  and  orientation  of  the  two  trajectories  are  presented  in  Table  5  for  all  ten 
runs. 


Figure  72:  Different  cell  trajectories  about  zero  point  and  relative  to  APG  local  grid.  Joint 
ArcSecond/IMU  data  in  red  and  IMU  only  data  in  black. 
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Figure  73:  Different  cell  trajectories  from  IMU  only  data  at  other  Arcsecond  measured  zero  points.  Local 

grid  coordinates. 
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Figure  74:  Comparison  of  IMU  (black)  and  ArcSecond  (red)  trajectories.  ArcSecond  coordinates. 
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Figure  75:  Second  comparison  of  IMU  (black)  and  ArcSecond  (red)  trajectories.  ArcSecond  coordinates. 
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Difference 
in  X  (m) 

Difference 
in  Y(m) 

Horizontal 

Difference 

(m) 

Difference 
in  Z  (m) 

Difference 
in  Yaw 
(degs) 

Difference 
in  Pitch 
(degs) 

Difference 
in  Roll 
(degs) 

Run  #1 

Mean 

-0.002 

0.001 

0.017 

-0.005 

1.96 

-2.41 

-1.37 

Standard  Deviation 

0.018 

0.013 

0.014 

0.006 

0.86 

0.61 

0.47 

Run  #2 

Mean 

-0.001 

-0.006 

0.013 

-0.003 

1.83 

-2.45 

-1.33 

Standard  Deviation 

0.011 

0.008 

0.009 

0.003 

0.74 

0.40 

0.35 

Run  #3 

Mean 

-0.004 

0.044 

0.046 

-0.004 

1.62 

-2.56 

-1.33 

Standard  Deviation 

0.007 

0.020 

0.017 

0.002 

0.62 

0.31 

0.27 

Run  #4 

Mean 

-0.002 

0.022 

0.027 

-0.003 

1.58 

-2.46 

-1.31 

Standard  Deviation 

0.012 

0.011 

0.007 

0.003 

0.60 

0.32 

0.31 

Run  #5 

Mean 

-0.003 

0.025 

0.028 

-0.002 

1.65 

-2.49 

-1.37 

Standard  Deviation 

0.007 

0.013 

0.009 

0.002 

0.73 

0.27 

0.28 

Run  #6 

Mean 

-0.003 

0.017 

0.020 

-0.003 

0.79 

-1.13 

-0.62 

Standard  Deviation 

0.007 

0.009 

0.008 

0.003 

0.64 

0.27 

0.39 

Run  #7 

Mean 

0.000 

0.000 

0.008 

-0.005 

0.85 

-1.04 

-0.66 

Standard  Deviation 

0.007 

0.007 

0.005 

0.002 

0.41 

0.19 

0.21 

Run  #8 

Mean 

-0.006 

0.026 

0.028 

-0.002 

0.73 

-1.21 

-0.71 

Standard  Deviation 

0.007 

0.011 

0.010 

0.005 

0.69 

0.56 

0.39 

Run  #9 

Mean 

0.000 

-0.001 

0.012 

-0.004 

0.94 

-1.09 

-0.62 

Standard  Deviation 

0.013 

0.010 

0.012 

0.003 

0.85 

0.38 

0.52 

Run  #10 

Mean 

0.011 

-0.024 

0.039 

0.005 

1.45 

-1.14 

-0.66 

Standard  Deviation 

0.027 

0.022 

0.021 

0.005 

0.67 

0.25 

0.29 

OVERALL 

Average  Mean 

-0.001 

0.010 

0.024 

-0.003 

1.341 

-1.798 

-0.997 

Average  Standard 
Deviation 

0.012 

0.012 

0.011 

0.003 

0.680 

0.356 

0.347 

Table  5:  Differences  between  ArcSecond  and  IMU  Trajectories. 


7.4.  Inversion  Results 

The  EM61-HH  data  was  collected  on  the  standard  Allegro  handheld  computer  provided  by 
Geonics.  The  data  was  collected  at  a  data  rate  of  15  Hz.  The  Geonics  handheld  has  a  second 
COM  port  to  collect  standardized  GPS  data.  ENSCO  is  working  on  using  this  port  to  provide 
time  stamping  between  the  IMU  data  and  the  EM61-HH  data.  However,  at  the  time  of  these  tests 
this  was  not  working  properly.  There  was  no  attempt  at  all  to  synchronize  the  ArcSecond  data  to 
the  Geonics  data.  The  EM61-HH  data  was  manually  shifted  to  within  a  second  of  the  other 
systems,  and  we  accounted  for  the  residual  time  difference  by  adding  it  in  as  a  fit  parameter  in 
our  EM  inversion  routine. 

The  basic  electromagnetic  induction  response  model  involves  a  simple  dipole  response  [6]. 
Inversion  of  this  model  solves  for  the  object’s  location,  object  orientation,  and  the  magnetic 
polarization  response  terms.  The  polarization  terms,  referred  to  here  as  P’s,  determine  the 
strength  of  the  object’s  induced  response  along  the  object’s  physical  axes.  It  is  these  terms  that 
allow  for  discrimination  between  buried  ordnance  and  other  metallic  clutter.  For  elongated 
ferrous  ordnance  items  such  as  a  projectile  or  mortar,  there  is  one  large,  primary  P  response 
along  the  item’s  long  axis  and  two  smaller,  secondary  responses  transverse  to  this,  Px  >  Py  =  pz- 
Accurately  determining  these  parameters  from  model-based  inversion  is  limited  by  the  signal-to- 
noise  ratio  (SNR)  of  the  sensor  data  and  by  the  positioning  accuracy  of  the  sensor  trajectory. 
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Simulations  performed  for  the  SERDP/ESTCP  Geolocation  Workshop  showed  that  positioning 
accuracy  of  better  than  1  centimeter  is  required  to  invert  P’s  with  less  than  20%  error. 

There  are  other  complications  in  inverting  the  EM61-HH  data  combined  with  the  IMU 
trajectory.  Over  minutes  of  data  collection,  the  zero  level  of  the  sensor  will  change.  The  data 
collected  here  starts  with  the  sensor  on  a  metal  tripod  with  an  unknown  offset.  It  then  moves 
back  and  forth  over  an  object  perhaps  reaching  zero  at  points,  perhaps  not.  Because  of  this,  an 
offset  parameter  has  been  added  to  the  data  inversion.  Also,  the  EM  data  is  not  synchronized 
with  the  IMU  data.  A  time  difference  tenn  is  yet  another  fit  parameter. 

Besides  the  fit  parameters,  the  inversion  algorithm  returns  a  reduced  chi-squared  value  for  the  fit 
and  a  squared  correlation  coefficient  R2  between  the  best  model  fit  and  the  measured  data.  The 
dipole  fit  error  is  expressed  as  V(l-f?2).  The  percentage  fit  error  grows  with  both  low  SNR 
signals  and  with  poor  positioning.  For  the  ten  separate  measurements  made  over  the  A3  cell,  the 
EM61-HH  data  was  inverted  first  using  the  ArcSecond  positioning  and  then  using  the  IMU 
positioning.  Figure  76  plots  the  fit  error  results  for  one  positioning  system  versus  the  other  using 
the  same  EM  data.  For  half  of  the  data,  the  two  systems  result  in  comparable  fit  errors.  For  the 
data  sets  where  the  trajectories  did  not  match  well,  the  IMU  fit  errors  are  significantly  greater. 
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Figure  76:  Fit  error  comparison  between  EM  data  with  IMU  and  ArcSecond  positioning. 


Figure  77  plots  the  spread  in  the  inverted  [1  parameters  over  all  20  mm  measurements.  In  these 
plots,  the  x  axis  plots  the  largest  [!  value  from  the  fit.  The  y  axis  plots  the  average  (diamond 
symbol)  and  the  range  (vertical  line)  of  the  two  secondary  P  values.  For  axisymmetric  objects 
like  the  20  mm,  the  secondary  P’s  should  be  equal  and  the  vertical  line  should  be  a  dot.  The 
points  in  77(a)  are  from  fits  using  Arcsecond  positioning.  The  mean  primary  P  value  is  0.387 
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with  a  standard  deviation  of  0.028  (7%  of  mean).  The  average  secondary  value  is  0.110  with  a 
standard  deviation  of  0.014  (13%).  The  same  fits  with  IMU  positioning  are  shown  in  12(b).  The 
mean  value  is  0.370  with  a  standard  deviation  of  0.059  (16%).  The  secondary  mean  is  0. 105  with 
a  standard  deviation  of  0.024  (23%).  The  spread  in  the  p  fits  is  roughly  twice  as  great  for  the 
IMU  positioning.  Similar  increased  spreads  in  the  other  object  fit  parameters  (X,  Y  location  and 
depth)  were  observed  with  the  IMU  positioning.  Figure  77(c)  plots  the  fit  results  for  IMU  only 
measurements  made  over  the  20  mm  cells:  A3,  C3,  and  D3.  The  results  are  comparable  to  77(b). 
Figure  78  plots  the  fit  results  for  data  collected  over  the  A4,  C4,  and  D4  cells  containing  a  40 
mm  grenade.  These  results  are  plotted  in  black  with  diamond  symbols.  For  comparison,  the  20 
mm  results  are  plotted  in  red.  There  are  comparable  spreads  in  this  data,  but  it  is  sufficiently 
small  to  discriminate  between  the  two  items.  It  should  be  noted  that  one  fit  result  is  significantly 
off  and  that  another  data  set  failed  to  fit  successfully. 

Figure  79  plots  the  fit  error  versus  the  SNR  of  the  EM61-FIH  signal  from  these  data  sets:  joint  20 
mm  ArcSecond  (green  X)  /  IMU  (red  X),  20  mm  IMU  only  (black  diamond),  and  40mm  grenade 
IMU  only  (blue  triangle).  The  general  trend  is  for  decreasing  fit  error  with  increasing  SNR.  The 
solid  black  curve  is  a  calculated  limit  for  fit  error  as  a  function  of  SNR  assuming  no  position 
errors.  The  dotted  curves  plot  the  limits  assuming  position  errors  on  the  order  of  one  and  two 
centimeters.  The  bulk  of  the  fits  fall  under  the  one  centimeter  curve. 

The  rest  of  the  data  collected  over  other  cells  on  the  APG  calibration  grid  had  SNR  in  the  range 
of  roughly  0  to  30  dB.  With  an  RMS  noise  level  of  ~4  millivolts,  this  is  in  the  range  of  a  few  to 
approximately  100  millivolts.  The  inversion  algorithm  had  trouble  converging  to  a  solution  for 
data  in  the  20  to  30  dB  range.  In  some  cases,  several  local  minima  were  found  by  the  fit 
algorithm. 

The  majority  of  the  joint  ArcSecond  and  ENSCO  IMU  positioning  data  showed  agreement  in  the 
horizontal  sensor  head  trajectories  to  be  within  one  to  three  centimeters  (see  Table  5).  There 
were  two  data  sets  where  the  ENSCO  IMU  data  was  found  to  drift  significantly  off  from  the 
ArcSecond  horizontal  positioning.  The  EM  fits  resulting  from  these  two  systems  showed  general 
agreement,  with  the  ArcSecond  showing  a  smaller  spread  in  the  fit  parameters  (7%  to  13% 
variations  about  the  mean  response  P’s  versus  16%  to  23%).  Even  with  the  larger  spread  in  P’s, 
the  IMU  positioned  fits  are  sufficient  to  discriminate  between  a  20  mm  projectile  and  a  40  mm 
grenade.  Comparisons  of  the  fit  results  overall  are  promising  with  the  bulk  of  the  fit  errors 
falling  about  theoretical  curves  indicating  roughly  one  centimeter  relative  positioning  error 
(figure  79). 

There  is  no  current  explanation  for  the  two  data  runs  that  drifted  significantly  off.  There  was 
nothing  indicated  in  the  IMU  output  to  flag  them  as  being  poor.  The  EM  fits  were  of  poorer,  but 
not  unreasonable,  quality  (10-14%  versus  5-10%).  The  fit  algorithm  did  have  trouble  converging 
for  this  data  and  there  were  several  P  parameter  outliers. 
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Figure  77:  Comparison  of  fitted  p  polarization  parameters  on  20  mm  cells. 
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Figure  78:  p  fits  of  40mm  grenade  data  compared  to  20  mm  fits. 


Figure  79:  Piots  of  fit  error  versus  EM61-HH  signai-to-noise  ratio.  Curves  indicate  fimits  based  on  0,  1, 
and  2  centimeter  positioning  errors. 
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8.  Adaptive  Processing  -  Joint  Target/Track  Estimation 


The  analysis  reported  in  the  previous  section  shows  that  while  the  Honeywell  IMU  is  capable  of 
providing  sub-centimeter  sensor  positioning  accuracy,  it  does  not  always  do  so.  As  part  of  this 
project  we  have  investigated  the  possibility  of  including  sensor  data  in  the  track  estimation 
process.  The  approach  is  based  on  the  principle  that  the  best  model/data  fit  occurs  when  the 
sensor  trajectory  is  correct.  Slight  errors  in  the  trajectory  degrade  the  fit  quality.  The  idea  is  then 
to  assume  that  the  trajectory  from  the  IMU  is  in  error,  and  allow  slight  modifications  of  the 
trajectory  in  the  EMI  inversion  process,  iterating  until  the  best  fit  to  the  data  is  achieved. 

Trajectory  errors  typically  arise  because  of  errors  in  estimating  the  biases  or  initial  and  final 
sensor  orientation.  From  §2  we  see  that  the  IMU  track  in  the  earth-fixed  frame  between  two  co¬ 
located  zero  velocity  updates  can  be  completely  specified  by  24  parameters.  If  5m,  5a ,  B  ,  and 
Ba  are  assumed  known  and  stable  (this  is  the  case  for  the  Honeywell  IMU),  we  are  left  with  12 

parameters  needed  to  specify  the  track  via  IMU  measurements  (with  the  implicit  assumption  that 
we  start  and  end  with  zero  velocity  updates  at  the  same  position).  These  are: 

A0)  -  the  initial  biases  in  the  angular  rate  measurements  (3); 

Aa  -  the  initial  biases  in  the  acceleration  measurements  (3); 

Si// , ,  S6j ,  S<j)i ,  Si// ,  ,  SO  f  and  S<f>  f  -  the  uncertainties  in  the  initial  and  final  yaw,  pitch  and  roll 
angles  (6) 

By  defining  a  function  that  accepts  these  12  parameters  as  an  input  vector  and  returns  a  measure 
of  goodness-of-fit,  we  can  use  a  downhill  simplex  method  (such  as  AMOEBA,  described  in 
section  10.4  of  Numerical  Recipes  in  C  [10])  to  perform  a  multidimensional  minimization  of  the 
function.  The  function  uses  the  12  input  parameters  to  construct  a  track  from  the  IMU  data,  thus 
associating  a  position  and  orientation  to  each  EM61-HH  reading,  and  uses  this  track  information 
as  input  data  our  standard  EM  Marquardt-Levenberg  inversion.  The  downhill  simplex  iterates  on 
track  parameter  hypotheses,  and  at  each  step,  the  EMI  inversion  gives  the  best  fit  target 
parameters.  The  correlation  squared,  R2,  is  computed  and  1-  R2  is  returned  as  the  goodness-of-fit 
measure  of  the  function. 

Figure  80  shows  results  of  applying  the  downhill  simplex  method  to  adaptively  refine  the  sensor 
trajectory  when  the  measured  trajectory  is  slightly  in  error.  It  is  based  on  simulated  EM  data  for 
an  object  with  the  following  characteristics: 

x0,  y0,  z0  =  (1.7,  2.7,  0.25)  m 

¥o,0o,(I>o=(21,2,  0)° 


Pi,  p2,  p3  =  (2,  0.5,  0.5) 


The  inset  diagram  shows  the  true  trajectory  in  black  and  the  input  (distorted)  trajectory  in  red. 
The  main  plot  in  Figure  80  shows  the  RMS  difference  between  the  true  trajectory  and  the 
estimated  trajectory  as  the  number  of  iterations  of  the  downhill  simplex  estimator  increases.  The 
deviation  of  the  initial  distorted  trajectory  from  the  true  trajectory  is  about  5  cm  rms.  The  search 
algorithm  converges  to  within  about  0.5  cm  rms  of  the  true  trajectory. 
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Figure  80:  Approach  to  true  trajectory  using  downhill  simplex  procedure.  Inset  shows  initial  distorted 
trajectory  in  red  and  true  trajectory  in  black. 


There  is  no  noise  in  the  simulation  shown  in  Figure  80.  We  have  repeated  the  procedure  with 
added  noise  of  2  mV  rms  and  4  mV  rms.  The  peak  signal  is  70  mV.  In  both  cases  the  procedure 
converges.  The  final  track  error  is  also  about  0.5  cm  rms  for  2  mV  noise  and  is  about  1  cm  rms 
for  4  mV  noise.  The  estimated  target  parameters  for  each  case  are  shown  in  Table  6.  There  is  a 
slight  degradation  in  performance  with  increasing  noise.  Figure  81  compares  the  observed 
variation  of  dipole  fit  error  with  noise  with  the  expected  variation.  The  dipole  fit  error  should  be 
equal  to  the  rms  noise  divided  by  the  rms  signal.  It  appears  that  there  is  a  slight  (few  %)  penalty 
for  adaptively  refining  the  track  estimate.  Note,  however,  that  with  the  original,  distorted  track 
the  dipole  fit  error  is  20%  in  the  noise-free  case. 


noise  (mV) 

Pi 

P2 

Ps 

depth  (m) 

fit  error  (%) 

ground  truth 

2.0 

0.50 

0.50 

0.250 

— 

0 

1.91 

0.48 

0.47 

0.249 

1.7 

2 

1.96 

0.53 

0.47 

0.252 

10.6 

4 

1.81 

0.51 

0.34 

0.245 

20.8 

Table  6:  Target  parameters  from  joint  target/track  estimation 
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Figure  82:  Fit  error  vs.  noise  level.  Dashed  line  corresponds  to  perfect  sensor  track  information. 
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9.  Conclusions 


The  objective  of  this  project  was  to  develop  an  inexpensive,  robust  way  to  determine  the  precise 
position  of  a  handheld  UXO  sensor  as  it  is  swept  about  a  suspected  buried  UXO  item,  enabling 
the  sensor  to  be  used  for  reliable  UXO/clutter  discrimination.  We  were  ale  to  demonstrate  that  a 
compact  (7  cm  square),  lightweight  (<400  gm)  Honeywell  HG1900  IMU  can  provide  sufficiently 
accurate  positioning  to  support  accurate  estimation  of  target  parameters  with  the  EM61-HH 
sensor.  The  Honeywell  unit  costs  about  $10K.  We  also  tested  a  less  expensive  (~$4K)  unit,  the 
Crossbow  IMU400CC-100,  but  found  that  it  did  not  provide  sufficiently  accurate  positioning. 
The  major  problem  with  the  Crossbow  unit  appears  to  be  that  its  angle  rate  noise  and  stability  are 
not  good  enough  to  accurately  track  attitude  and  keep  gravitational  acceleration  from 
contaminating  the  horizontal  acceleration  channels.  The  angular  rate  random  walk  spec  for  the 
Crossbow  unit  is  2.25  deg/hr  ,  while  the  corresponding  spec  for  the  Honeywell  unit  is  0.1 
deg/hr1/2. 

During  the  course  of  the  project,  we  found  that  the  response  characteristics  of  the  EM61-HH  can 
significantly  affect  the  accuracy  of  target  parameter  estimates  when  the  sensor  is  moving  about 
over  the  target.  This  occurs  because  the  sensor  output  is  integrated  with  an  analog  filter  that  both 
shifts  and  distorts  the  sensor's  response.  We  measured  the  filter  response  characteristics  and 
included  it  in  a  dynamic  forward  model  for  the  EM61-HH.  In  order  to  get  discrimination-quality 
accuracy  we  have  to  use  the  dynamic  model  when  inverting  data  collected  as  the  sensor  is  swept 
about  over  a  target. 

In  a  few  of  the  test  cases  with  the  Honeywell  we  had  an  uncorrected  drift  of  several  cm  relative 
to  ground  truth  provided  by  a  laser  (ArcSecond)  positioning  system.  We  are  able  correct  for  this 
sort  of  drift  using  an  adaptive  processing  approach  which  includes  the  EMI  sensor  data  in  the 
track  estimation  process.  The  trajectory  errors  arise  because  of  errors  in  estimating  the  IMU 
biases.  By  including  the  bias  parameters  in  our  EM  inversion  procedure,  we  can  simultaneously 
determine  the  best  track  estimate  along  with  the  target  parameters  (location  and  principal  axis 
polarizabilities).  Simulations  using  this  procedure  show  that  we  can  reduce  5  cm  rms  track  errors 
to  <5  mm  rms. 
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Appendix  A:  The  Yaw-Pitch-Roll  Rotation  Sequence 


The  rotation  matrix  R bo^b(0  that  transforms  the  initial  body  frame  at  time  t  =  to  to  the  body 

frame  at  time  t  is  defined  by  an  ordered  sequence  of  three  plane  rotations  involving  the  Euler 
angles  i//,  # and  (f>  (yaw,  pitch  and  roll,  respectively). 


YAW  (\p)  -  rotation  about  PITCH  (0)  -  rotation  about  ROLL  (<|))  -  rotation  about 

the  z-axis  the  new  y-axis  y'  the  last  x-axis  x" 

Figure  Al:  The  yaw-pitch-roll  rotation  sequence. 


The  yaw  rotation  aligns  the  new  x'  -axis  with  the  projection  of  the  body  x-axis  at  time  t  into  the 
initial  body  frame. 


(Al) 

f 

X 

r 

y 

= 

cost//  sin^  0 
-sin^  cos^  0 

X 

y 

=R, 

X 

y 

z' 

0  0  1 

z 

z 

The  pitch  rotation  aligns  the  new  x"  -axis  with  the  body  x-axis  at  time  t. 


(A2) 

rr 

X 

y" 

= 

cos  6  0  -  sin  6 

0  1  0 

r 

X 

y' 

=  Re 
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X 

y' 

rr 

Z 

sin#  0  cos# 

r 

Z 

r 

Z 

The  roll  rotation  aligns  the  new  y"'  -axis  and  z"  -axis  with  the  body  y-axis  and  z-axis, 
respectively,  at  time  t. 
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Thus,  vectors  represented  in  the  initial  body-frame  at  time  t  =  to  are  transformed  into  the  body- 
frame  at  time  t  by  the  sequence  of  three  rotations: 

(A4) 

where 


Ra0->*  = 


cos  y/  cos  # 

-  sin  y/  cos  (j)  +  cos  y/  sin  #  sin  (f> 
sin  y/  sin  (j)  +  cos  y/  sin  #  cos  (j) 


sin  y/  cos  6  -  sin  6 

cos  y/  cos  ^  +  sin  y/  sin  9  sin  (j)  cos  9  sin  (j) 
-  cos  y/  sin  (j)  +  sin  y/  sin  9  cos  (j)  cos  9  cos  (j) 


The  inverse  transformation  is  (recalling  the  fact  that  because  R  is 

orthonormal) 


(A5) 


v  °  =  R,  ,  v  =  R 


i.e. 


R 


b^>b0 


cos^cos  9 
sin  cos# 
-  sin# 


-  sin  y/  cos  (/)  +  cos  y/  sin  9  sin  (j> 
cos  y/  cos  (j)  +  sin  y/  sin  9  sin  (j> 
cos#  sin  ^ 


sin  y/  sin  (f>  +  cos  y/  sin  #  cos  (j) 
-  cos  y/  sin  (j>  +  sin  y/  sin  #  cos  (j) 
cos#  cos  ^ 
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